text:
"""
All RadeonProRender types inherit from RPRObject
"""
abstract type RPRObject{PtrType} end

Base.pointer(object::RPRObject) = object.pointer
Base.cconvert(::Type{T}, object::RPRObject{T}) where {T<:Ptr} = object
Base.unsafe_convert(::Type{T}, object::RPRObject{T}) where {T} = pointer(object)
function Base.unsafe_convert(::Type{Ptr{Nothing}}, object::RPRObject{T}) where {T}
    return Base.unsafe_convert(Ptr{Nothing}, pointer(object))
end

function release(object::RPRObject)
    remove_from_context!(object.parent, object)
    if object.pointer != C_NULL
        rprObjectDelete(object)
    end
    object.pointer = C_NULL
    if hasproperty(object, :referenced_memory)
        object.referenced_memory = nothing
    end
    return
end

@enum PluginType Tahoe Northstar Hybrid HybridPro

using Libdl

function plugin_path(plugin::PluginType)
    if plugin == Tahoe
        return RadeonProRender_jll.libTahoe64
    elseif plugin == Northstar
        return RadeonProRender_jll.libNorthstar64
    end
    Sys.isapple() && error("Only Tahoe and Northstar are supported")
    path_base = dirname(RadeonProRender_jll.libNorthstar64)
    # Hybrid isn't currently part of Binarybuilder products, since its missing on Apple
    if plugin == Hybrid
        return joinpath(path_base, "Hybrid.$(Libdl.dlext)")
    elseif plugin == HybridPro
        return joinpath(path_base, "HybridPro.$(Libdl.dlext)")
    end
end

mutable struct Context <: RPRObject{rpr_context}
    pointer::rpr_context
    objects::Base.IdSet{RPRObject}
    plugin::PluginType
    function Context(plugin::PluginType, creation_flags, singleton=true, props=C_NULL, cache_path=C_NULL)
        id = RPR.rprRegisterPlugin(plugin_path(plugin))
        @assert(id != -1)
        plugin_ids = [id]
        ctx_ptr = RPR.rprCreateContext(RPR.RPR_API_VERSION, plugin_ids, 1, creation_flags, props, cache_path)
        RPR.rprContextSetActivePlugin(ctx_ptr, id)
        ctx = new(ctx_ptr, Base.IdSet{RPRObject}(), plugin)
        if singleton
            if isassigned(ACTIVE_CONTEXT)
                @info("releasing old context")
                release(ACTIVE_CONTEXT[])
            end
            ACTIVE_CONTEXT[] = ctx
        end
        return ctx
    end
end

const ACTIVE_CONTEXT = Base.RefValue{Context}()

get_context(object::RPRObject) = get_context(object.parent)
get_context(context::Context) = context

function add_to_context!(contextlike::RPRObject, object::RPRObject)
    return push!(get_context(contextlike).objects, object)
end

function remove_from_context!(contextlike::RPRObject, object::RPRObject)
    return delete!(get_context(contextlike).objects, object)
end

"""
    Context()
Empty constructor, defaults to using opencl, GPU0 and Tahoe as the rendering plugin.

* `resource::rpr_creation_flags_t = RPR_CREATION_FLAGS_ENABLE_GPU0`
* `plugin::PluginType=Northstar`: can be Tahoe, Northstar or Hybrid or HybridPro. Everything is tested with Northstar, which is a complete rewrite of Tahoe. Hybrid is for real time rendering and only supports Uber material.
* `singleton = true`:  set to true, to only allow one context at the same time. This is useful, to immediately free all old resources when creating a new context, instead of relying on the GC.
"""
function Context(; plugin=Northstar, resource=RPR_CREATION_FLAGS_ENABLE_GPU0, singleton=true)
    return Context(plugin, resource, singleton)
end

"""
Wraps a RadeonProRender type into a gc registered higher level type.
arguments:
`name`: name of the Julia type
`typ`: type of the RadeonProRender object
`target`: RadeonProRender constructor function
`args`: arguments to constructor function
`super` Julia super type
"""
macro rpr_wrapper_type(name, typ, target, args, super)
    argnames = args.args
    parent = argnames[1]
    expr = quote
        mutable struct $name{PT} <: $super{$typ}
            pointer::$typ
            parent::PT
            $name(pointer::$typ, parent::PT) where {PT} = new{PT}(pointer, parent)
            function $name($(argnames...))
                ptr = $target($(argnames...))
                object = new{typeof($parent)}(ptr, $parent)
                add_to_context!($parent, object)
                return object
            end
        end
    end
    return esc(expr)
end

function release(context::Context)
    if context.pointer == C_NULL
        @assert isempty(context.objects)
        return
    end
    scenes = Scene[]
    for object in copy(context.objects)
        if object isa Scene
            push!(scenes, object)
        else
            release(object)
        end
    end
    empty!(context.objects)
    foreach(s -> (rprSceneClear(s); release(s)), scenes)
    empty!(scenes)
    rprContextClearMemory(context)
    rprObjectDelete(context)
    context.pointer = C_NULL
    return
end

@rpr_wrapper_type MaterialSystem rpr_material_system rprContextCreateMaterialSystem (context, typ) RPRObject
@rpr_wrapper_type MaterialNode rpr_material_node rprMaterialSystemCreateNode (matsystem, typ) RPRObject
@rpr_wrapper_type Camera rpr_camera rprContextCreateCamera (context,) RPRObject
@rpr_wrapper_type Scene rpr_scene rprContextCreateScene (context,) RPRObject
@rpr_wrapper_type PostEffect rpr_post_effect rprContextCreatePostEffect (context, type) RPRObject

function set!(pe::PostEffect, param::String, x::Number)
    return rprPostEffectSetParameter1f(pe, param, x)
end

mutable struct VoxelGrid <: RPRObject{rpr_grid}
    pointer::rpr_grid
    parent::Context
    referenced_memory::Any
end

function VoxelGrid(context::Context, nx, ny, nz, grid_values::AbstractArray, grid_indices::AbstractArray)
    grid_ptr = Ref{rpr_grid}()

    grid_indices = convert(Vector{UInt64}, vec(grid_indices))
    valuesf0 = convert(Vector{Float32}, vec(grid_values))

    rprContextCreateGrid(context, grid_ptr, nx, ny, nz, grid_indices, length(grid_indices),
                         RPR.RPR_GRID_INDICES_TOPOLOGY_I_U64, valuesf0, sizeof(valuesf0), 0)
    grid = VoxelGrid(grid_ptr[], context, (valuesf0, grid_indices))
    add_to_context!(context, grid)
    return grid
end

function VoxelGrid(context, values::AbstractArray{<:Number,3})
    indices = UInt64.((1:length(values)) .- 1)
    return VoxelGrid(context, size(values)..., values, indices)
end

mutable struct Curve <: RPRObject{rpr_curve}
    pointer::rpr_curve
    parent::Context
    referenced_memory::Any
end

function Curve(context::Context, control_points::AbstractVector, indices::AbstractVector,
               radius::AbstractVector, uvs::AbstractVector, segments_per_curve::AbstractVector,
               creationFlag_tapered=0)
    curve_ref = Ref{rpr_curve}()
    control_pointsf0 = convert(Vector{Point3f}, control_points)
    indices_i = convert(Vector{RPR.rpr_int}, indices)
    segments_per_curve_i = convert(Vector{RPR.rpr_int}, segments_per_curve)
    uvs_f0 = convert(Vector{Vec2f}, uvs)
    radiusf0 = convert(Vector{Float32}, radius)
    rprContextCreateCurve(context, curve_ref, length(control_pointsf0), control_pointsf0, sizeof(Point3f),
                          length(indices_i), length(segments_per_curve_i), indices_i, radiusf0, uvs_f0,
                          segments_per_curve_i, creationFlag_tapered)
    curve = Curve(curve_ref[], context, (control_pointsf0, indices_i, segments_per_curve_i, uvs_f0, radiusf0))
    add_to_context!(context, curve)
    return curve
end

@rpr_wrapper_type HeteroVolume rpr_hetero_volume rprContextCreateHeteroVolume (context,) RPRObject

function set_albedo_lookup!(volume::HeteroVolume, colors::AbstractVector)
    return RPR.rprHeteroVolumeSetAlbedoLookup(volume, convert(Vector{RGB{Float32}}, colors), length(colors))
end

function set_albedo_grid!(volume::HeteroVolume, albedo::VoxelGrid)
    return RPR.rprHeteroVolumeSetAlbedoGrid(volume, albedo)
end

function set_density_lookup!(volume::HeteroVolume, density::AbstractVector)
    return RPR.rprHeteroVolumeSetDensityLookup(volume, convert(Vector{Vec3f}, density), length(density))
end

function set_density_grid!(volume::HeteroVolume, density::VoxelGrid)
    return RPR.rprHeteroVolumeSetDensityGrid(volume, density)
end

function Base.push!(scene::Scene, volume::HeteroVolume)
    return rprSceneAttachHeteroVolume(scene, volume)
end

mutable struct Shape <: RPRObject{rpr_shape}
    pointer::rpr_shape
    parent::Context
    referenced_arrays::Any
end

function set!(shape::Shape, volume::HeteroVolume)
    return rprShapeSetHeteroVolume(shape, volume)
end

"""
Abstract AbstractLight Type
"""
abstract type AbstractLight{PtrType} <: RPRObject{PtrType} end

@rpr_wrapper_type EnvironmentLight rpr_light rprContextCreateEnvironmentLight (context,) AbstractLight
@rpr_wrapper_type PointLight rpr_light rprContextCreatePointLight (context,) AbstractLight
@rpr_wrapper_type DirectionalLight rpr_light rprContextCreateDirectionalLight (context,) AbstractLight
@rpr_wrapper_type SkyLight rpr_light rprContextCreateSkyLight (context,) AbstractLight
@rpr_wrapper_type SpotLight rpr_light rprContextCreateSpotLight (context,) AbstractLight

"""
Default shape constructor which works with every Geometry from the package
GeometryTypes (Meshes and geometry primitives alike).
"""
function Shape(context::Context, meshlike; kw...)
    m = uv_normal_mesh(meshlike; kw...)
    v, n, fs, uv = decompose(Point3f, m), normals(m), faces(m), texturecoordinates(m)
    if isnothing(n)
        n = normals(v, fs)
    end
    return Shape(context, v, n, fs, uv)
end

#=
The yield + nospecialize somehow prevent some stack/memory corruption when run with `--check-bounds=yes`
This is kind of magical and still under investigation.
MWE:
https://gist.github.com/SimonDanisch/475064ae102141554f65e926f3070630
=#
@nospecialize
function Shape(context::Context, vertices, normals, faces, uvs)
    @assert length(vertices) == length(normals)
    @assert length(vertices) == length(uvs)

    vraw = decompose(Point3f, vertices)
    nraw = decompose(Vec3f, normals)
    uvraw = map(uv -> Vec2f(1 - uv[2], 1 - uv[1]), uvs)
    f = decompose(TriangleFace{OffsetInteger{-1,rpr_int}}, faces)
    iraw = collect(reinterpret(rpr_int, f))
    facelens = fill(rpr_int(3), length(faces))

    @assert eltype(vraw) == Point3f
    @assert eltype(nraw) == Vec3f
    @assert eltype(uvraw) == Vec2f

    foreach(i -> checkbounds(vertices, i + 1), iraw)
    yield() # grrr, why you!!
    rpr_mesh = rprContextCreateMesh(context, vraw, length(vertices), sizeof(Point3f), nraw, length(normals),
                                    sizeof(Vec3f), uvraw, length(uvs), sizeof(Vec2f), iraw, sizeof(rpr_int),
                                    iraw, sizeof(rpr_int), iraw, sizeof(rpr_int), facelens, length(faces))

    jl_references = (vraw, nraw, uvraw, iraw, facelens)
    shape = Shape(rpr_mesh, context, jl_references)
    add_to_context!(context, shape)
    return shape
end
@specialize

"""
Creating a shape from a shape is interpreted as creating an instance.
"""
function Shape(context::Context, shape::Shape)
    inst = Shape(rprContextCreateInstance(context, shape), context, nothing)
    add_to_context!(context, inst)
    return inst
end

"""
Create layered shader give two shaders and respective IORs
"""
function layeredshader(matsys, base, top, weight=(0.5, 0.5, 0.5, 1.0))
    layered = MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_BLEND)

    set!(layered, RPR.RPR_MATERIAL_INPUT_COLOR0, base)
    # Set shader for top layer
    set!(layered, RPR.RPR_MATERIAL_INPUT_COLOR1, top)
    # Set index of refraction for base layer
    set!(layered, RPR.RPR_MATERIAL_INPUT_WEIGHT, weight...)
    return layered
end

"""
FrameBuffer wrapper type
"""
mutable struct FrameBuffer <: RPRObject{rpr_framebuffer}
    pointer::rpr_framebuffer
    parent::Context
end

function FrameBuffer(context::Context, c::Type{C}, dims::NTuple{2,Int}) where {C<:Colorant}
    desc = Ref(RPR.rpr_framebuffer_desc(dims...))
    fmt = RPR.rpr_framebuffer_format(length(c), RPR_COMPONENT_TYPE_FLOAT32)
    frame_buffer = FrameBuffer(RPR.rprContextCreateFrameBuffer(context, fmt, desc), context)
    add_to_context!(context, frame_buffer)
    return frame_buffer
end

function get_data(framebuffer::FrameBuffer)
    framebuffer_size = Ref{Cint}()
    RPR.rprFrameBufferGetInfo(framebuffer, RPR.RPR_FRAMEBUFFER_DATA, 0, C_NULL, framebuffer_size)
    data = zeros(RGBA{Float32}, framebuffer_size[] ÷ sizeof(RGBA{Float32}))
    @assert sizeof(data) == framebuffer_size[]
    RPR.rprFrameBufferGetInfo(framebuffer, RPR.RPR_FRAMEBUFFER_DATA, framebuffer_size[], data, C_NULL)
    return data
end

"""
Image wrapper type
"""
mutable struct Image <: RPRObject{rpr_image}
    pointer::rpr_image
    parent::Context
    referenced_memory::Any
end

"""
Automatically generated the correct `rpr_image_format` for an array
"""
rpr_image_format(a::Array) = rpr_image_format(eltype(a))

function rpr_image_format(::Type{T}) where {T<:Union{Colorant,Number}}
    return RPR.rpr_image_format(T <: Number ? 1 : length(T), component_type(T))
end

"""
Gets the correct component type for Images from array element types
"""
component_type(x::Type{T}) where {T<:Colorant} = component_type(eltype(T))
component_type(x::Type{T}) where {T<:Union{UInt8,Colors.N0f8}} = RPR_COMPONENT_TYPE_UINT8
component_type(x::Type{T}) where {T<:Union{Float16}} = RPR_COMPONENT_TYPE_FLOAT16
component_type(x::Type{T}) where {T<:Union{Float32}} = RPR_COMPONENT_TYPE_FLOAT32

"""
Creates the correct `rpr_image_desc` for a given array.
"""
function rpr_image_desc(image::Array{T,N}) where {T,N}
    row_pitch = size(image, 1) * sizeof(T)
    return RPR.rpr_image_desc(ntuple(i -> N < i ? 0 : size(image, i), 3)..., row_pitch, 0)
end

"""
Automatically creates an Image from a matrix of colors
"""
function Image(context::Context, image::Array{T,N}) where {T,N}
    desc_ref = Ref(rpr_image_desc(image))
    img = Image(rprContextCreateImage(context, rpr_image_format(image), desc_ref, image), context,
                (image, desc_ref))
    add_to_context!(context, img)
    return img
end

function Image(context::Context, image::AbstractArray)
    return Image(context, collect(image))
end

"""
Automatically loads an image from the given `path`
"""
function Image(context::Context, path::AbstractString)
    img = Image(rprContextCreateImageFromFile(context, path), context, nothing)
    add_to_context!(context, img)
    return img
end

"""
Sets the radiant power for a given Light
"""
function setradiantpower!(light::AbstractLight, color::Colorant{T,3}) where {T}
    c = RGB{Float32}(color)
    return setradiantpower!(light, comp1(c), comp2(c), comp3(c))
end

function setradiantpower!(light::PointLight, r::Number, g::Number, b::Number)
    return rprPointLightSetRadiantPower3f(light, rpr_float(r), rpr_float(g), rpr_float(b))
end

function setradiantpower!(light::SpotLight, r::Number, g::Number, b::Number)
    return rprSpotLightSetRadiantPower3f(light, rpr_float(r), rpr_float(g), rpr_float(b))
end

function setradiantpower!(light::DirectionalLight, r::Number, g::Number, b::Number)
    return rprDirectionalLightSetRadiantPower3f(light, rpr_float(r), rpr_float(g), rpr_float(b))
end

"""
Sets the image for an EnvironmentLight
"""
function set!(light::EnvironmentLight, image::Image)
    return rprEnvironmentLightSetImage(light, image)
end

function set!(light::EnvironmentLight, context::Context, image::Array)
    return set!(light, Image(context, image))
end

function set!(light::EnvironmentLight, context::Context, path::AbstractString)
    return set!(light, Image(context, path))
end

"""
Sets the intensity scale an EnvironmentLight
"""
function setintensityscale!(light::EnvironmentLight, intensity_scale::Number)
    return rprEnvironmentLightSetIntensityScale(light, intensity_scale)
end

"""
Sets a portal for a Light
"""
function setportal!(light::AbstractLight, portal::AbstractGeometry)
    return setportal!(light, Shape(portal))
end

function setportal!(light::EnvironmentLight, portal::Shape)
    return rprEnvironmentLightSetPortal(light, portal)
end

function setportal!(light::SkyLight, portal::Shape)
    return rprSkyLightSetPortal(light, portal)
end

"""
Sets the albedo for a SkyLight
"""
setalbedo!(skylight::SkyLight, albedo::Number) = rprSkyLightSetAlbedo(skylight.x, albedo)

"""
Sets the turbidity for a SkyLight
"""
setturbidity!(skylight::SkyLight, turbidity::Number) = rprSkyLightSetTurbidity(skylight.x, turbidity)

"""
Sets the scale for a SkyLight
"""
setscale!(skylight::SkyLight, scale::Number) = rprSkyLightSetScale(skylight, scale)

function set!(context::Context, parameter::rpr_context_info, f::AbstractFloat)
    return rprContextSetParameterByKey1f(context, parameter, f)
end

function set!(context::Context, parameter::rpr_context_info, ui::Integer)
    return rprContextSetParameterByKey1u(context, parameter, ui)
end
function set!(context::Context, parameter::rpr_context_info, ui)
    # overload for enums
    return rprContextSetParameterByKey1u(context, parameter, ui)
end

set!(context::Context, aov::rpr_aov, fb::FrameBuffer) = rprContextSetAOV(context, aov, fb)

function set!(shape::Shape, image::Image)
    return rprShapeSetDisplacementImage(shape, image)
end

function set!(shape::Shape, context::Context, image::Array)
    return set!(shape, Image(context, image))
end

"""
Sets arbitrary values to RadeonProRender objects
"""
set!(context::Context, framebuffer::FrameBuffer) = rprContextSetFrameBuffer(context, framebuffer)

function set!(base::MaterialNode, parameter::rpr_material_node_input, a::Number, b::Number, c::Number,
              d::Number)
    return rprMaterialNodeSetInputFByKey(base, parameter, a, b, c, d)
end

function set!(base::MaterialNode, parameter::rpr_material_node_input, ui::Integer)
    return rprMaterialNodeSetInputUByKey(base, parameter, ui)
end

function set!(base::MaterialNode, parameter::rpr_material_node_input, f::Vec4)
    return set!(base, parameter, f...)
end

function set!(base::MaterialNode, parameter::rpr_material_node_input, color::Colorant)
    c = RGBA{Float32}(color)
    return set!(base, parameter, comp1(c), comp2(c), comp3(c), alpha(c))
end

function set!(mat::MaterialNode, parameter::rpr_material_node_input, enum)
    return set!(mat, parameter, UInt(enum))
end

function set!(shape::Shape, material::MaterialNode)
    return rprShapeSetMaterial(shape, material)
end

function set!(shape::Curve, material::MaterialNode)
    return rprCurveSetMaterial(shape, material)
end

function set!(material::MaterialNode, parameter::rpr_material_node_input, material2::MaterialNode)
    return rprMaterialNodeSetInputNByKey(material, parameter, material2)
end

function set!(material::MaterialNode, parameter::rpr_material_node_input, image::Image)
    return rprMaterialNodeSetInputImageDataByKey(material, parameter, image)
end

function set!(material::MaterialNode, parameter::rpr_material_node_input, context::Context, image::Array)
    return set!(material, parameter, Image(context, image))
end

function set!(context::Context, scene::Scene)
    return rprContextSetScene(context, scene)
end

function set!(scene::Scene, light::EnvironmentLight)
    return rprSceneSetEnvironmentLight(scene, light)
end

function set!(scene::Scene, camera::Camera)
    return rprSceneSetCamera(scene, camera)
end

"""
Sets the displacement scale for a shape
"""
function setdisplacementscale!(shape::Shape, scale::Number)
    return rprShapeSetDisplacementScale(shape, rpr_float(scale))
end

"""
Sets the subdivions for a shape
"""
function setsubdivisions!(shape::Shape, subdivions::Integer)
    return rprShapeSetSubdivisionFactor(shape, rpr_uint(subdivions))
end

"""
Transforms firerender objects with a transformation matrix
"""
function transform!(shape::Shape, transform::AbstractMatrix)
    return rprShapeSetTransform(shape, false, convert(Matrix{Float32}, transform))
end

function transform!(light::AbstractLight, transform::AbstractMatrix)
    return rprLightSetTransform(light, false, convert(Matrix{Float32}, transform))
end

function transform!(camera::Camera, transform::AbstractMatrix)
    return rprCameraSetTransform(camera, false, convert(Matrix{Float32}, transform))
end

"""
Clears a scene
"""
clear!(scene::Scene) = rprSceneClear(scene)

"""
Clears a framebuffer
"""
clear!(frame_buffer::FrameBuffer) = rprFrameBufferClear(frame_buffer)

"""
Pushes objects to Scene
"""
function Base.push!(scene::Scene, shape::Shape)
    return rprSceneAttachShape(scene, shape)
end

function Base.push!(scene::Scene, curve::Curve)
    return rprSceneAttachCurve(scene, curve)
end

function Base.push!(scene::Scene, light::AbstractLight)
    return rprSceneAttachLight(scene, light)
end

function set!(context::Context, pe::PostEffect)
    return rprContextAttachPostEffect(context, pe)
end

function Base.delete!(scene::Scene, light::AbstractLight)
    return rprSceneDetachLight(scene, light)
end

function Base.delete!(scene::Scene, shape::Shape)
    return rprSceneDetachShape(scene, shape)
end

"""
Gets the currently attached EnvironmentLight from a scene
"""
function Base.get(scene::Scene, ::Type{EnvironmentLight})
    return EnvironmentLight(rprSceneGetEnvironmentLight(scene))
end

"""
Gets the currently attached Camera from a scene
"""
function Base.get(scene::Scene, ::Type{Camera})
    return Camera(rprSceneGetCamera(scene))
end

"""
Sets the lookat of a camera
"""
function lookat!(camera::Camera, position::Vec3, lookatvec::Vec3, upvec::Vec3)
    return rprCameraLookAt(camera, position..., lookatvec..., upvec...)
end

"""
renders the context
"""
function render(context::Context)
    return rprContextRender(context)
end

"""
Renders a tile of the context
"""
function render(context::Context, tile::Rect)
    mini, maxi = minimum(tile), maximum(tile)
    return rprContextRenderTile(context, mini[1], maxi[1], mini[2], maxi[2])
end

"""
Saves a picture of the current framebuffer at `path`
"""
function save(fb::FrameBuffer, path::AbstractString)
    return rprFrameBufferSaveToFile(fb, path)
end

#=
Some iterators for tiled rendering
=#

"""
Iterates randomly over an iterable
"""
struct RandomIterator{T}
    iterable::T
end

function Base.iterate(ti::RandomIterator, state=rand(1:length(ti.iterable)))
    return (ti.iterable[state], rand(1:length(ti.iterable)))
end

"""
Iterates randomly over all elements in an iterable.
The middle has a higher sampling probability.
"""
struct MiddleSampleIterator{T}
    iterable::T
end

function multiplier(i, j, n1, n2)
    a = abs(mod(i - div(n1, 2), n1) - div(n1, 2))
    b = abs(mod(j - div(n2, 2), n2) - div(n2, 2))
    return max(a, b)
end

function init_sample_pool(ti)
    samplepool = Tuple{Int,Int}[]
    for i in 1:size(ti.iterable, 1)
        for j in 1:size(ti.iterable, 2)
            append!(samplepool, fill((i, j), multiplier(i, j, size(ti.iterable)...)))
        end
    end
    return samplepool
end

function Base.iterate(ti::MiddleSampleIterator, state=init_sample_pool(ti))
    isempty(state) && return nothing
    i = rand(1:length(state))
    i2 = state[i]
    filter!(state) do s
        return s != i2
    end
    i = sub2ind(size(ti.iterable), i2...)
    return (ti.iterable[i], state)
end

"""
Iterates over all tiles in a Rectangular parent tile.
"""
struct TileIterator
    size::Vec{2,Int}
    tile_size::Vec{2,Int}
    lengths::Vec{2,Int}
end

function TileIterator(size, tile_size)
    s, ts = Vec(size), Vec(tile_size)
    lengths = Vec{2,Int}(ceil(Vec{2,Float64}(s) ./ Vec{2,Float64}(ts)))
    return TileIterator(s, ts, lengths)
end

function Base.getindex(ti::TileIterator, i, j)
    ts = Vec(ti.tile_size)
    xymin = (Vec(i, j) .- 1) .* ts
    xymax = xymin .+ ts
    return Rect(xymin, min(xymax, ti.size) .- xymin)
end

function Base.getindex(ti::TileIterator, linear_index::Int)
    i, j = ind2sub(size(ti), linear_index)
    return ti[i, j]
end

Base.size(ti::TileIterator, i) = ti.lengths[i]
Base.size(ti::TileIterator) = (ti.lengths...,)
Base.length(ti::TileIterator) = prod(ti.lengths)

function Base.iterate(ti::TileIterator, state=1)
    length(ti) > state && return nothing
    return ti[state], state + 1
end

"""
Customizable defaults for the most common tonemapping operations
"""
function set_standard_tonemapping!(context; typ=RPR.RPR_TONEMAPPING_OPERATOR_PHOTOLINEAR,
                                   photolinear_sensitivity=0.5f0, photolinear_exposure=0.5f0,
                                   photolinear_fstop=2.0f0, reinhard02_prescale=1.0f0,
                                   reinhard02_postscale=1.0f0, reinhard02_burn=1.0f0, linear_scale=1.0f0,
                                   aacellsize=4.0, imagefilter_type=RPR.RPR_FILTER_BLACKMANHARRIS,
                                   aasamples=4.0)
    norm = PostEffect(context, RPR.RPR_POST_EFFECT_NORMALIZATION)
    set!(context, norm)
    tonemap = PostEffect(context, RPR.RPR_POST_EFFECT_TONE_MAP)
    set!(context, tonemap)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_TYPE, typ)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_LINEAR_SCALE, linear_scale)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_PHOTO_LINEAR_SENSITIVITY, photolinear_sensitivity)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_PHOTO_LINEAR_EXPOSURE, photolinear_exposure)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_PHOTO_LINEAR_FSTOP, photolinear_fstop)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_REINHARD02_PRE_SCALE, reinhard02_prescale)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_REINHARD02_POST_SCALE, reinhard02_postscale)
    set!(context, RPR.RPR_CONTEXT_TONE_MAPPING_REINHARD02_BURN, reinhard02_burn)
    set!(context, RPR.RPR_CONTEXT_IMAGE_FILTER_TYPE, imagefilter_type)

    gamma = PostEffect(context, RPR.RPR_POST_EFFECT_GAMMA_CORRECTION)
    set!(context, gamma)
    set!(context, RPR.RPR_CONTEXT_DISPLAY_GAMMA, 1.0)
    # set!(context, "aacellsize", aacellsize)
    # set!(context, "aasamples", aasamples)
    # set!(context, RPR.RPR_CONTEXT_AA_ENABLED, UInt(1))

    # bloom = PostEffect(context, RPR.RPR_POST_EFFECT_BLOOM)
    # set!(context, bloom)
    # set!(bloom, "weight", 1.0f0)
    # set!(bloom, "radius", 0.4f0)
    # set!(bloom, "threshold", 0.2f0)
    return
end


text:
# Math functions unrelated to any particular structure.

@inline lerp(t::AbstractFloat, v₁::AbstractFloat, v₂::AbstractFloat) =
    (1.0 - t) * v₁ + t * v₂

@inline Γ(n::T) where T<:AbstractFloat = (n * eps(T)) / (1 - n * eps(T))

ext: jl


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module CarmichaelNumbers

using Nemo, NumberTheory, Counts

export ModuleCarmichaelNumbers
export isCarmichael, I002997, F002997, L002997
export isweakCarmichael, I225498, F225498, L225498

"""
* isCarmichael, I002997, F002997, L002997
* isweakCarmichael, I225498, F225498, L225498
"""
const ModuleCarmichaelNumbers = ""

"""
Is ``n`` a Carmichael/Šimerka number?
"""
function isCarmichael(n)
    (n == 1 || isEven(n) || isPrime(n)) && return false
    for f in Factors(n)
        (f[2] > 1 || (n - 1) % (f[1] - 1) != 0) && return false
    end
    return true
end

"""
Iterate over the first n Carmichael/Šimerka numbers.
"""
I002997(n) = takeFirst(isCarmichael, n)

"""
Iterate over the Carmichael/Šimerka numbers which do not exceed n.
"""
F002997(n) = filter(isCarmichael, 1:n)

"""
Return the first n Carmichael/Šimerka numbers in an array.
"""
L002997(n) = collect(I002997(n))

"""
Is ``n`` a weak Carmichael number?
"""
function isweakCarmichael(n)
    (n == 1 || isEven(n) || isPrime(n)) && return false
    for f in Factors(n)
        (n - 1) % (f[1] - 1) != 0 && return false
    end
    return true
end

"""
Iterate over the first n weak Carmichael numbers.
"""
I225498(n) = takeFirst(isweakCarmichael, n)

"""
Iterate over the weak Carmichael numbers which do not exceed n.
"""
F225498(n) = filter(isweakCarmichael, 1:n)

"""
Return the first n weak Carmichael numbers in an array.
"""
L225498(n) = collect(I225498(n))

#START-TEST-########################################################

using Test, SeqTests

function test()
    @testset "Carmichael" begin
        @test ! isCarmichael(560)
        @test isCarmichael(561)
        @test ! isCarmichael(563)
        @test isweakCarmichael(561)
        @test !isweakCarmichael(563)
        @test isweakCarmichael(625)

        if is_oeis_installed()
            L = [L002997, L225498]
            SeqTest(L, 'L')
        end

    end
end

function demo()
    for n in 1:30000
        isCarmichael(n) && println(n)
    end

    println()
    for v in I225498(10)
        println(v)
    end
end

"""
@time L = L002997(30)
0.967192 seconds (4.01 M allocations: 248.863 MiB, 28.97% gc time)
"""
function perf()
    @time L = L002997(30)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
# SPDX-License-Identifier: MIT

#======================================================================
This file wrapped a few basic linear algebra operations in order to
have more fine-grained control between CPU/GPU code and different
semifield algebra operations.
======================================================================#

## Element-wise division.

eldiv!(out, x::AbstractArray, y::AbstractArray) = broadcast!(/, out, x, y)

## Element-wise multiplication.

elmul!(out, x::AbstractArray, y::AbstractArray) =  broadcast!(*, out, x, y)
function elmul!(out, x::AbstractSparseArray{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    fill!(out, zero(T))
    broadcast!(*, out, x, y)
end
elmul!(out, x::CuSparseVector, y::AbstractVector) = elmul_svdv!(out, x, y)
elmul!(out, x::CuSparseVector, Y::AbstractMatrix) = elmul_svdm!(out, x, Y)

## Matrix multiplication.

function matmul!(out, X::AbstractMatrix{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    mul!(out, X, y, one(T), zero(T))
end

function matmul!(out, X::AbstractSparseMatrix{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    fill!(out, zero(T))
    mul!(out, X, y, one(T), zero(T))
end

matmul!(out, X::CuMatrix, Y::CuMatrix) = mul_dmdm!(out, X, Y)
matmul!(out, X::CuSparseMatrixCSR, y::AbstractVector) = mul_smdv!(out, X, y)
matmul!(out, X::CuSparseMatrixCSR, Y::AbstractMatrix) = mul_smdm!(out, X, Y)


text:
function Base.summary(io::IO, A::AbstractDimArray{T,N}) where {T,N}
    if N == 0
        print(io, "0-dimensional ")
    elseif N == 1
        print(io, size(A, 1), "-element ")
    else
        print(io, join(size(A), "×"), " ")
    end
    print(io, string(nameof(typeof(A)), "{$T,$N}"))
end

function Base.show(io::IO, mime::MIME"text/plain", A::AbstractDimArray{T,N}) where {T,N}
    lines = 0
    summary(io, A)
    print_name(io, name(A))
    lines += Dimensions.print_dims(io, mime, dims(A))
    !(isempty(dims(A)) || isempty(refdims(A))) && println(io)
    lines += Dimensions.print_refdims(io, mime, refdims(A))
    println(io)

    # Printing the array data is optional, subtypes can
    # show other things here instead.
    ds = displaysize(io)
    ioctx = IOContext(io, :displaysize => (ds[1] - lines, ds[2]))
    show_after(ioctx, mime, A)

    return nothing
end

# Semi-interface methods for adding addional `show` text
# for AbstractDimArray/AbstractDimStack subtypes
# TODO actually document in the interface
show_after(io::IO, mime, A::AbstractDimArray) = print_array(io, mime, parent(A))

# Showing the array is optional for AbstractDimArray
# `print_array` must be called from `show_after`.
function print_array(io::IO, mime, A::AbstractArray{T,0}) where T
    print(_print_array_ctx(io, T), "\n", A[])
end
function print_array(io::IO, mime, A::AbstractArray{T,1}) where T
    Base.print_matrix(_print_array_ctx(io, T), A)
end
function print_array(io::IO, mime, A::AbstractArray{T,2}) where T
    Base.print_matrix(_print_array_ctx(io, T), A)
end
function print_array(io::IO, mime, A::AbstractArray{T,N}) where {T,N}
    o = ones(Int, N-2)
    frame = A[:, :, o...]
    onestring = join(o, ", ")
    println(io, "[:, :, $(onestring)]")
    Base.print_matrix(_print_array_ctx(io, T), frame)
    nremaining = prod(size(A, d) for d=3:N) - 1
    nremaining > 0 && print(io, "\n[and ", nremaining," more slices...]")
end

function _print_array_ctx(io, T)
    IOContext(io, :compact=>true, :limit=>true, :typeinfo=>T)
end
# print a name of something, in yellow
function print_name(io::IO, name)
    if !(name == Symbol("") || name isa NoName)
        printstyled(io, string(" ", name); color=:yellow)
    end
end


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadde32 = QuadDE(Float32)
quadde64 = QuadDE(Float64)
quaddeBF = QuadDE(BigFloat)

# Test integrals are cited from:
# Kahaner, D.K.: Comparison of numerical quadrature formulas, Mathematical
# software, Rice, J.R. (Ed.), Academic Press. pp.229-259 (1971)

# NOTE that since I couldn't get the original book, the problems below are
# cited from a second source.
# http://id.nii.ac.jp/1001/00011109/

# Test problem 1
let
    f(x::AbstractFloat) = exp(x)
    expect = BigFloat("1.71828182845904523536028747135")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 2
# non-smooth function
let
    f(x::AbstractFloat) = floor(min(x/3*10, one(x)))
    expect = BigFloat("7.00000000000000000000000000000e-1")

    # Unfortunately, this problem is very difficult for the default QuadDE.
    # The accuracy is (slowly) depending only on the step size of trapezoidal rule.

    # In practical use, user may not know what actually f is.
    # However, if you know where, splitting the integral interval at the
    # discontinuity may help.
    I, E = quadde32(f, 0, 0.3, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 0.3, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-18
    I, E = quaddeBF(f, 0, 0.3, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 3
let
    f(x::AbstractFloat) = sqrt(x)
    expect = BigFloat("6.66666666666666666666666666666e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 4
let
    f(x::T) where {T<:AbstractFloat} = cosh(x)*92/100 - cos(x)
    expect = BigFloat("4.7942822668880166735857796183531e-1")

    I, E = quadde32(f, -1, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -1, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 5
let
    f(x::T) where {T<:AbstractFloat} = 1/(x^4 + x^2 + T(0.9))
    f(x::BigFloat) = 1/(x^4 + x^2 + BigFloat("0.9"))
    expect = BigFloat("1.5822329637296729331174689490262e0")

    I, E = quadde32(f, -1, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -1, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 6
let
    f(x::AbstractFloat) = x*sqrt(x)
    expect = BigFloat("4.00000000000000000000000000000e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 7
# Singular point at x = 0
let
    f(x::AbstractFloat) = 1/sqrt(x)
    expect = 2

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 8
let
    f(x::AbstractFloat) = 1/(x^4 + 1)
    expect = BigFloat("8.66972987339911037573995163882e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 9
let
    f(x::AbstractFloat) = 2/(2 + sin(x*314159/10000))
    expect = BigFloat("1.15470066904371304340692220986e0")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 10
let
    f(x::AbstractFloat) = 1/(1 + x)
    expect = BigFloat("6.93147180559945309417232121458e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 11
let
    f(x::AbstractFloat) = 1/(exp(x) + 1)
    expect = BigFloat("3.7988549304172247536823662649e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 12
# Singular point at x = 0
let
    f(x::AbstractFloat) = x/(exp(x) - 1)
    expect = BigFloat("7.77504634112248276417586545425e-1")

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 13
# Attenuating oscillation
let
    f(x::AbstractFloat) = sin(x*314159/1000)/(x*314159/100000)
    expect = BigFloat("9.0986452565692970698e-3")

    I, E = quadde32(f, 0.1, 1)
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0.1, 1)
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-18
    I, E = quaddeBF(f, 0.1, 1, rtol=rtol)
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 14
# Super fast decay around x ~ 0
# FIXME: split the integral range into parts
let
    f(x::T) where {T<:AbstractFloat} = 5*sqrt(T(2))*exp(x^2*(-50)*314_159/100_000)
    expect = BigFloat("5.0000021117e-1")

    I, E = quadde32(f, 0, 0.5, 10)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 0.5, 10)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-10
    I, E = quaddeBF(f, 0, 0.5, 10, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol*2)
    @test E ≤ rtol*norm(I)
end

# Test problem 15
# Super fast decay around x ~ 0
# FIXME: split the integral range into parts
let
    f(x::AbstractFloat) = 25*exp(-x*25)
    expect = 1

    I, E = quadde32(f, 0, 1, 10)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1, 10)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, 10, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol*2)
    @test E ≤ rtol*norm(I)
end

# Test problem 16
# Super fast decay around x ~ 0
let
    f(x::AbstractFloat) = 1/(2500*x^2 + 1)*50/314_159*100_000
    expect = BigFloat("4.99363802871016550828171090341e-1")

    I, E = quadde32(f, 0, 10)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 10)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 10, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 17
# Attenuating oscillation
let
    f(x::AbstractFloat) = (sin(x*50*314_159/100_000)/(x*50*314_159/100_000))^2*50
    expect = BigFloat("1.12139569626709460839885e-1")

    I, E = quadde32(f, 0.01, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0.01, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    # NOTE: the accuracy is not improved so much?
    rtol = 1e-16
    I, E = quaddeBF(f, 0.01, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 18
# Oscillating function
# NOTE: expect value is doubtful?
let
    f(x::AbstractFloat) = cos(cos(x) + 3*sin(x) + 2*cos(2x) + 3*sin(2x) + 3*cos(3x))
    expect = BigFloat("8.386763233809718250439e-1")

    I, E = quadde32(f, 0, π)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, π)
    @test I isa Float64
    @test_broken I ≈ expect
    @test_skip E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-11
    I, E = quaddeBF(f, 0, π, rtol=rtol)
    @test I isa BigFloat
    @test_broken isapprox(I, expect, rtol=rtol)
    @test_skip E ≤ rtol*norm(I)
end

# Test problem 19
let
    f(x::AbstractFloat) = log(x)
    expect = -1

    I, E = quadde32(f, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 20
let
    f(x::T) where {T<:AbstractFloat} = 1/(x^2 + T(1.005))
    f(x::BigFloat) = 1/(x^2 + BigFloat("1.005"))
    expect = BigFloat("1.5643964440690497730914930158085e0")

    I, E = quadde32(f, -1, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -1, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test problem 21
# including spikes
# FIXME: split the integral range into parts
let
    f(x::AbstractFloat) = 1/cosh((10x-2))^2+1/cosh((100x-40))^4+1/cosh((1000x-600))^6
    expect = BigFloat("0.2108027355005492773756")

    I, E = quadde32(f, 0, 0.4, 0.6, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 0.4, 0.6, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    atol = 1e-22
    I, E = quaddeBF(f, 0, 0.4, 0.6, 1, atol=atol)
    @test I isa BigFloat
    @test isapprox(I, expect, atol=atol)
    @test E ≤ atol
end


text:
module ReactiveMPRuleTest

using Test
using ReactiveMP
using MacroTools

import ReactiveMP: rule_macro_parse_on_tag, rule_macro_parse_fn_args, call_rule_macro_parse_fn_args

@testset "rule" begin

    @testset "rule_macro_parse_on_tag(expression)" begin
        @test rule_macro_parse_on_tag(:(:out)) == (:(Type{ Val{ :out } }), nothing, nothing)
        @test rule_macro_parse_on_tag(:(:mean)) == (:(Type{ Val{ :mean } }), nothing, nothing)
        @test rule_macro_parse_on_tag(:(:mean, k)) == (:(Tuple{ Val{ :mean }, Int }), :k, :(k = on[2]))
        @test rule_macro_parse_on_tag(:(:precision, r)) == (:(Tuple{ Val{ :precision }, Int }), :r, :(r = on[2]))

        @test_throws ErrorException rule_macro_parse_on_tag(:(out))
        @test_throws ErrorException rule_macro_parse_on_tag(:(123))
        @test_throws ErrorException rule_macro_parse_on_tag(:(:mean, 1))
        @test_throws ErrorException rule_macro_parse_on_tag(:(precision, r))
    end

    @testset "rule_macro_parse_fn_args(inputs; specname, prefix, proxy)" begin
        names, types, init = rule_macro_parse_fn_args([ (:m_out, :PointMass), (:m_mean, :NormalMeanPrecision) ]; specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))

        @test names == :(Type{ Val{ (:out, :mean) } })
        @test types == :(Tuple{ ReactiveMP.Message{ <: PointMass }, ReactiveMP.Message{ <: NormalMeanPrecision } })
        @test init == Expr[:(m_out = getdata(messages[1])), :(m_mean = getdata(messages[2]))]

        names, types, init = rule_macro_parse_fn_args([ (:m_out, :PointMass), (:m_mean, :NormalMeanPrecision) ]; specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

        @test names == :Nothing
        @test types == :Nothing
        @test init == Expr[]

        names, types, init = rule_macro_parse_fn_args([ (:m_out, :PointMass), (:q_mean, :NormalMeanPrecision) ]; specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

        @test names == :(Type{ Val{ (:mean, ) } })
        @test types == :(Tuple{ ReactiveMP.Marginal{ <: NormalMeanPrecision }, })
        @test init == Expr[ :(q_mean = getdata(marginals[1])) ]
    end

    @testset "call_rule_macro_parse_fn_args(inputs; specname, prefix, proxy)" begin
        names, values = call_rule_macro_parse_fn_args([ (:m_out, :(PointMass(1.0))), (:m_mean, :(NormalMeanPrecision(0.0, 1.0))) ]; specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))

        @test names == :(Val{ (:out, :mean) })
        @test values == :(ReactiveMP.Message(PointMass(1.0), false, false), ReactiveMP.Message(NormalMeanPrecision(0.0, 1.0), false, false))

        names, values = call_rule_macro_parse_fn_args([ (:m_out, :(PointMass(1.0))), (:m_mean, :(NormalMeanPrecision(0.0, 1.0))) ]; specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

        @test names == :nothing
        @test values == :nothing

        names, values = call_rule_macro_parse_fn_args([ (:m_out, :(PointMass(1.0))), (:q_mean, :(NormalMeanPrecision(0.0, 1.0))) ]; specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

        @test names == :(Val{ (:mean, ) })
        @test values == :((ReactiveMP.Marginal(NormalMeanPrecision(0.0, 1.0), false, false),))
    end

end

end


text:
# X and Z components in fT
# conducting 100 Ωm
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   10, 10, 10]

ross_primary = [25.7847     -18.7313]
ross_secondary = [7.7025     -8.6221
  5.9817     -7.5021
   5.049      -6.804
  4.1531     -6.0555
  3.2195     -5.1767
    2.36     -4.2461
  1.6135     -3.3041
  1.0446     -2.4585
0.64716     -1.7593
0.38013      -1.202
0.20975    -0.77834
0.10887    -0.47774
0.05337    -0.27878
0.02495     -0.1561
0.010711   -0.081757]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
import Plots, Plots.plot, Plots.annotate!, Plots.current, Plots.text, Plots.font

function plot(grid::TreeGrid{2}; markers::Bool = false, max_marker_level::Int = 5, path::Bool = false, indices = false, filter::Function = active)
    FullyThreadedTree.plot(grid.tree, markers = markers, max_marker_level = max_marker_level, path = path, filter = filter)
    if indices
        max_level = levels(grid.tree) - 1
        font_size = 10.0/(1<<min(max_level, max_marker_level))
        for cell ∈ cells(grid.tree, filter = filter, max_level = max_marker_level)
            pos = centroid(cell)
            annotate!(pos[1], pos[2], index(cell), font(font_size))
        end
        for face ∈ FullyThreadedTree.all_faces(grid.tree, filter = filter, max_level = max_marker_level)
            pos = centroid(face)
            annotate!(pos[1], pos[2], index(face), font(font_size))
        end
    end
    Plots.current()
end


text:
include("layers.jl")
include("cost_functions.jl")
include("nn.jl")

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
xavail = trues(15,15)
xcnt = 0
xgoal = 10
println("Testing adder with 4-bit words:")
while xcnt < xgoal
    m = rand(1:15)
    n = rand(1:15)
    xavail[m,n] || continue
    xavail[m,n] = xavail[n,m] = false
    xcnt += 1
    (s, c) = adder(m, n)
    oflow = c ? "*" : ""
    print(@sprintf "    %2d + %2d = %2d => " m n m+n)
    println(@sprintf("%s + %s = %s%s",
                     bits(m)[end-3:end],
                     bits(n)[end-3:end],
                     bits(s), oflow))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Driver script for all unit tests

using Jive

runtests(@__DIR__, skip=["runtests.jl", "helper", "broken", "_manual_"])


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SimplificationMethod

A method for simplifying geometric objects.
"""
abstract type SimplificationMethod end

"""
    simplify(object, method)

Simplify `object` with given `method`.
"""
function simplify end

function simplify(box::Box{2}, method::SimplificationMethod)
  c = simplify(boundary(box), method)
  PolyArea(c)
end

function simplify(polygon::Polygon, method::SimplificationMethod)
  c = [simplify(chain, method) for chain in chains(polygon)]
  PolyArea(c[1], c[2:end])
end

function simplify(multi::Multi, method::SimplificationMethod)
  Multi([simplify(geom, method) for geom in multi])
end

function simplify(domain::Domain, method::SimplificationMethod)
  GeometrySet([simplify(elem, method) for elem in domain])
end

"""
    decimate(object, ϵ)

Simplify `object` with an appropriate
simplification method and tolerance `ϵ`.
"""
decimate(object, ϵ) = simplify(object, DouglasPeucker(ϵ))

# ----------------
# IMPLEMENTATIONS
# ----------------

include("simplification/douglaspeucker.jl")


text:
"""
    Elkan()

Elkan algorithm implementation, based on "Charles Elkan. 2003.
Using the triangle inequality to accelerate k-means.
In Proceedings of the Twentieth International Conference on
International Conference on Machine Learning (ICML’03). AAAI Press, 147–153."

This algorithm provides much faster convergence than Lloyd algorithm especially
for high dimensional data.
It can be used directly in `kmeans` function

```julia
X = rand(30, 100_000)   # 100_000 random points in 30 dimensions

kmeans(Elkan(), X, 3) # 3 clusters, Elkan algorithm
```
"""
struct Elkan <: AbstractKMeansAlg end

function kmeans!(alg::Elkan, containers, X, k;
                n_threads = Threads.nthreads(),
                k_init = "k-means++", max_iters = 300,
                tol = 1e-6, verbose = false, init = nothing)
    nrow, ncol = size(X)
    centroids = init == nothing ? smart_init(X, k, n_threads, init=k_init).centroids : deepcopy(init)

    update_containers(alg, containers, centroids, n_threads)
    @parallelize n_threads ncol chunk_initialize(alg, containers, centroids, X)

    converged = false
    niters = 0
    J_previous = 0.0

    # Update centroids & labels with closest members until convergence
    while niters < max_iters
        niters += 1
        # Core iteration
        @parallelize n_threads ncol chunk_update_centroids(alg, containers, centroids, X)

        # Collect distributed containers (such as centroids_new, centroids_cnt)
        # in paper it is step 4
        collect_containers(alg, containers, n_threads)

        J = sum(containers.ub)

        # auxiliary calculation, in paper it's d(c, m(c))
        calculate_centroids_movement(alg, containers, centroids)

        # lower and ounds update, in paper it's steps 5 and 6
        @parallelize n_threads ncol chunk_update_bounds(alg, containers, centroids)

        # Step 7, final assignment of new centroids
        centroids .= containers.centroids_new[end]

        if verbose
            # Show progress and terminate if J stopped decreasing.
            println("Iteration $niters: Jclust = $J")
        end

        # Check for convergence
        if (niters > 1) & (abs(J - J_previous) < (tol * J))
            converged = true
            break
        end

        # Step 1 in original paper, calulation of distance d(c, c')
        update_containers(alg, containers, centroids, n_threads)
        J_previous = J
    end

    @parallelize n_threads ncol sum_of_squares(containers, X, containers.labels, centroids)
    totalcost = sum(containers.sum_of_squares)

    # Terminate algorithm with the assumption that K-means has converged
    if verbose & converged
        println("Successfully terminated with convergence.")
    end

    # TODO empty placeholder vectors should be calculated
    # TODO Float64 type definitions is too restrictive, should be relaxed
    # especially during GPU related development
    return KmeansResult(centroids, containers.labels, Float64[], Int[], Float64[], totalcost, niters, converged)
end

function create_containers(::Elkan, k, nrow, ncol, n_threads)
    lng = n_threads + 1
    centroids_new = Vector{Array{Float64,2}}(undef, lng)
    centroids_cnt = Vector{Vector{Int}}(undef, lng)

    for i = 1:lng
        centroids_new[i] = zeros(nrow, k)
        centroids_cnt[i] = zeros(k)
    end

    centroids_dist = Matrix{Float64}(undef, k, k)

    # lower bounds
    lb = Matrix{Float64}(undef, k, ncol)

    # upper bounds
    ub = Vector{Float64}(undef, ncol)

    # r(x) in original paper, shows whether point distance should be updated
    stale = ones(Bool, ncol)

    # distance that centroid moved
    p = Vector{Float64}(undef, k)

    labels = zeros(Int, ncol)

    # total_sum_calculation
    sum_of_squares = Vector{Float64}(undef, n_threads)

    return (
        centroids_new = centroids_new,
        centroids_cnt = centroids_cnt,
        labels = labels,
        centroids_dist = centroids_dist,
        lb = lb,
        ub = ub,
        stale = stale,
        p = p,
        sum_of_squares = sum_of_squares
    )
end

function chunk_initialize(::Elkan, containers, centroids, X, r, idx)
    ub = containers.ub
    lb = containers.lb
    centroids_dist = containers.centroids_dist
    labels = containers.labels
    centroids_new = containers.centroids_new[idx]
    centroids_cnt = containers.centroids_cnt[idx]

    @inbounds for i in r
        min_dist = distance(X, centroids, i, 1)
        label = 1
        lb[label, i] = min_dist
        for j in 2:size(centroids, 2)
            # triangular inequality
            if centroids_dist[j, label] > min_dist
                lb[j, i] = min_dist
            else
                dist = distance(X, centroids, i, j)
                label = dist < min_dist ? j : label
                min_dist = dist < min_dist ? dist : min_dist
                lb[j, i] = dist
            end
        end
        ub[i] = min_dist
        labels[i] = label
        centroids_cnt[label] += 1
        for j in axes(X, 1)
            centroids_new[j, label] += X[j, i]
        end
    end
end

function update_containers(::Elkan, containers, centroids, n_threads)
    # unpack containers for easier manipulations
    centroids_dist = containers.centroids_dist

    k = size(centroids_dist, 1) # number of clusters
    @inbounds for j in axes(centroids_dist, 2)
        min_dist = Inf
        for i in j + 1:k
            d = distance(centroids, centroids, i, j)
            centroids_dist[i, j] = d
            centroids_dist[j, i] = d
            min_dist = min_dist < d ? min_dist : d
        end
        for i in 1:j - 1
            min_dist = min_dist < centroids_dist[j, i] ? min_dist : centroids_dist[j, i]
        end
        centroids_dist[j, j] = min_dist
    end

    # TODO: oh, one should be careful here. inequality holds for eucledian metrics
    # not square eucledian. So, for Lp norm it should be something like
    # centroids_dist = 0.5^p. Should check one more time original paper
    centroids_dist .*= 0.25

    return centroids_dist
end

function chunk_update_centroids(::Elkan, containers, centroids, X, r, idx)
    # unpack
    ub = containers.ub
    lb = containers.lb
    centroids_dist = containers.centroids_dist
    labels = containers.labels
    stale = containers.stale
    centroids_new = containers.centroids_new[idx]
    centroids_cnt = containers.centroids_cnt[idx]

    @inbounds for i in r
        label_old = labels[i]
        label = label_old
        min_dist = ub[i]
        # tighten the loop, exclude points that very close to center
        min_dist <= centroids_dist[label, label] && continue
        for j in axes(centroids, 2)
            # tighten the loop once more, exclude far away centers
            j == label && continue
            min_dist <= lb[j, i] && continue
            min_dist <= centroids_dist[j, label] && continue

            # one calculation per iteration is enough
            if stale[i]
                min_dist = distance(X, centroids, i, label)
                lb[label, i] = min_dist
                ub[i] = min_dist
                stale[i] = false
            end

            if (min_dist > lb[j, i]) | (min_dist > centroids_dist[j, label])
                dist = distance(X, centroids, i, j)
                lb[j, i] = dist
                if dist < min_dist
                    min_dist = dist
                    label = j
                end
            end
        end

        if label != label_old
            labels[i] = label
            centroids_cnt[label_old] -= 1
            centroids_cnt[label] += 1
            for j in axes(X, 1)
                centroids_new[j, label_old] -= X[j, i]
                centroids_new[j, label] += X[j, i]
            end
        end
    end
end

function collect_containers(alg::Elkan, containers, n_threads)
    if n_threads == 1
        @inbounds containers.centroids_new[end] .= containers.centroids_new[1] ./ containers.centroids_cnt[1]'
    else
        @inbounds containers.centroids_new[end] .= containers.centroids_new[1]
        @inbounds containers.centroids_cnt[end] .= containers.centroids_cnt[1]
        @inbounds for i in 2:n_threads
            containers.centroids_new[end] .+= containers.centroids_new[i]
            containers.centroids_cnt[end] .+= containers.centroids_cnt[i]
        end

        @inbounds containers.centroids_new[end] .= containers.centroids_new[end] ./ containers.centroids_cnt[end]'
    end
end

function calculate_centroids_movement(alg::Elkan, containers, centroids)
    p = containers.p
    centroids_new = containers.centroids_new[end]

    for i in axes(centroids, 2)
        p[i] = distance(centroids, centroids_new, i, i)
    end
end

function chunk_update_bounds(alg, containers, centroids, r, idx)
    p = containers.p
    lb = containers.lb
    ub = containers.ub
    stale = containers.stale
    labels = containers.labels

    @inbounds for i in r
        for j in axes(centroids, 2)
            lb[j, i] = lb[j, i] > p[j] ? lb[j, i] + p[j] - 2*sqrt(abs(lb[j, i]*p[j])) : 0.0
        end
        stale[i] = true
        ub[i] += p[labels[i]] + 2*sqrt(abs(ub[i]*p[labels[i]]))
    end
end