text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module CSG

using FromFile: @from
using StaticArrays: SVector, SMatrix
using LinearAlgebra: norm

@from "../constants.jl" import VIS_RECURSION_LIMIT,
    MIN_VIS_TRI_SIZE, MESH_PRECISION
@from "Surfaces.jl" import Surfaces:
    Primitive,
    Surface,
    surfaceintersection,
    normal,
    interface,
    ParametricSurface,
    point,
    partials,
    uvrange,
    inside,
    onsurface,
    uv,
    samplesurface,
    TriangleMesh,
    Triangle,
    vertex,
    validtri,
    makemesh

@from "OpticalInterfaces.jl" import OpticalInterfaces: AllOpticalInterfaces
@from "BoundingBoxes.jl" import BoundingBoxes:
    BoundingBox, doesintersect, intersection, BoundingBoxes
@from "Transforms.jl" import Transforms: Transform, identitytransform
@from "Rays.jl" import Rays: AbstractRay, Ray, α
@from "Intersections.jl" import Intersections:
    Intersection,
    EmptyInterval,
    Interval,
    DisjointUnion,
    closestintersection,
    intervalintersection,
    RayOrigin,
    Infinity,
    threadedintervalpool,
    IntervalPool,
    emptyintervalpool!,
    reversenormal,
    intervalcomplement,
    lower,
    upper,
    intervalunion

# export CSGTree,
#     CSGGenerator,
#     leaf,
#     Intersection,
#     Infinity,
#     RayOrigin,
#     EmptyInterval,
#     Interval,
#     DisjointUnion,
#     closestintersection

# ============================================================================ #
#                             dP""b8 .dP"Y8  dP""b8
#                            dP   `" `Ybo." dP   `"
#                            Yb      o.`Y8b Yb  "88
#                             YboodP 8bodP'  YboodP
# ============================================================================ #

"""Abstract type representing any evaluated CSG structure."""
abstract type CSGTree{T} <: Primitive{T} end

"""
    ComplementNode{T,C<:CSGTree{T}} <: CSGTree{T}

An evaluated complement node within the CSG tree, must be the second child of a
[`IntersectionNode`](@ref) forming a subtraction.
"""
struct ComplementNode{T, C <: CSGTree{T}} <: CSGTree{T}
    child::C

    function ComplementNode(child::C) where {T <: Real, C <: CSGTree{T}}
        return new{T, C}(child)
    end
end
function Base.show(io::IO, a::ComplementNode{T}) where {T}
    return print(io, "Complement($(a.child))")
end
BoundingBox(a::ComplementNode{T}) where {T <: Real} = BoundingBox(a.child)

"""
    UnionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}

An evaluated union node within the CSG tree.
"""
struct UnionNode{T, L <: CSGTree{T}, R <: CSGTree{T}} <: CSGTree{T}
    leftchild::L
    rightchild::R
    bbox::BoundingBox{T}

    function UnionNode(
        a::L,
        b::R,
    ) where {T <: Real, L <: CSGTree{T}, R <: CSGTree{T}}
        # union should never contain a complement so should be fine
        # FIXME it errors out here
        return new{T, L, R}(a, b, union(BoundingBox(a), BoundingBox(b)))
    end
end
function Base.show(io::IO, a::UnionNode{T}) where {T}
    return print(io, "Union($(a.leftchild), $(a.rightchild))")
end

"""
    IntersectionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}

An evaluated intersection node within the CSG tree.
"""
struct IntersectionNode{T, L <: CSGTree{T}, R <: CSGTree{T}} <: CSGTree{T}
    leftchild::L
    rightchild::R
    bbox::BoundingBox{T}

    function IntersectionNode(
        a::L,
        b::R,
    ) where {T <: Real, L <: CSGTree{T}, R <: CSGTree{T}}
        # normal intersection is fine for most nodes
        return new{T, L, R}(a, b, intersection(BoundingBox(a), BoundingBox(b)))
    end

    function IntersectionNode(
        a::L,
        b::R,
    ) where {T <: Real, L <: CSGTree{T}, R <: ComplementNode{T}}
        # this is a subtraction so just take the original unclipped bounding box for simplicity
        return new{T, L, R}(a, b, BoundingBox(a))
    end
end
function Base.show(io::IO, a::IntersectionNode{T}) where {T}
    return print(io, "Intersection($(a.leftchild), $(a.rightchild))")
end

"""
    LeafNode{T,S<:ParametricSurface{T}} <: CSGTree{T}

An evaluated leaf node in the CSG tree, `geometry` attribute which contains a
[`ParametricSurface`](@ref) of type `S`.
The leaf node also has a transform associated which is the composition of all
nodes above it in the tree.
As such, transforming points from the geometry using this transform puts them
in world space, and transforming rays by the inverse transform puts them in
object space.
"""
struct LeafNode{T, S <: ParametricSurface{T, 3}} <: CSGTree{T}
    geometry::S
    transform::Transform{T}
    invtransform::Transform{T}
    bbox::BoundingBox{T}

    function LeafNode(
        a::S,
        transform::Transform{T},
    ) where {T <: Real, S <: ParametricSurface{T, 3}}
        # store the transformed bounding box so nodes higher in the tree have correct global space bounding boxes
        return new{T, S}(
            a,
            transform,
            inv(transform),
            BoundingBox(a, transform),
        )
    end
end
function Base.show(io::IO, a::LeafNode{T}) where {T}
    if a.transform == identitytransform(T)
        print(io, "Leaf($(a.geometry))")
    else
        print(io, "Leaf($(a.geometry), $(a.transform))")
    end
end

BoundingBox(a::CSGTree{T}) where {T <: Real} = a.bbox

"""
    CSGGenerator{T<:Real}

This is the type you should use when making CSG objects.
This type allows for the construction of [`CSGTree`](@ref) objects with
different transforms.
When the generator is evaluated, all transforms are propagated down to the
[`LeafNode`](@ref)s and stored there.

# Example
```julia
a = Cylinder(1.0,1.0)
b = Plane([0.0,0.0,1.0], [0.0,0.0,0.0])
generator = a ∩ b
# now make a csg object that can be ray traced
csgobj = generator(Transform(1.0,1.0,2.0))
```
"""
struct CSGGenerator{T <: Real}
    f::Function
end

function (a::CSGGenerator{T})(
    transform::Transform{T},
)::CSGTree{T} where {T <: Real}
    return a.f(transform)
end
function (a::CSGGenerator{T})()::CSGTree{T} where {T <: Real}
    return a.f(identitytransform(T))
end

"""
    leaf(surf::ParametricSurface{T}, transform::Transform{T} =
        identitytransform(T)) -> CSGGenerator{T}

Create a leaf node from a parametric surface with a given transform.
"""
function leaf(
    surf::ParametricSurface{T},
    transform::Transform{T} = identitytransform(T),
) where {T <: Real}
    return CSGGenerator{T}(
        (parenttransform) -> LeafNode(surf, parenttransform * transform),
    )
end

"""
    transform(surf::CSGGenerator{T}, transform::Transform{T} =
        identitytransform(T)) -> CSGGenerator{T}

Returns a new CSGGenerator with another transform applied. This is useful if
you want multiple copies of a premade CSG structure with different transforms,
for example in an MLA.
"""
function transform(
    n::CSGGenerator{T},
    transform::Transform{T} = identitytransform(T),
) where {T <: Real}
    return CSGGenerator{T}((parenttransform) -> n(parenttransform * transform))
end

# CSG objects are trees, not graphs, since it makes no sense to reuse
# intermediate results. Hence no need to keep track of common subexpressions.
# Static arrays are much faster than mutable arrays so want node transforms to
# be static. Unfortunately the transformations cascade from the root to the
# leaves but the nodes have to be created from the leaves to the roots.
# Wrap the csg operations in function that delays evaluation until the transform
# has been computed.

"""
    ∩(a::Union{CSGGenerator{T},ParametricSurface{T}}, b::Union{CSGGenerator{T},
        ParametricSurface{T}}) where {T<:Real}

Create a binary node in the CSG tree representing an intersection between `a` and `b`.

![Intersect Image](https://upload.wikimedia.org/wikipedia/commons/0/0b/Boolean_intersect.PNG)
"""
Base.:∩

function Base.:∩(a::CSGGenerator{T}, b::CSGGenerator{T}) where {T <: Real}
    return CSGGenerator{T}(
        (parenttransform) ->
            IntersectionNode(a(parenttransform), b(parenttransform)),
    )
end
Base.:∩(a::ParametricSurface, b::ParametricSurface) = leaf(a) ∩ leaf(b)
Base.:∩(a::ParametricSurface, b::CSGGenerator) = leaf(a) ∩ b
Base.:∩(a::CSGGenerator, b::ParametricSurface) = a ∩ leaf(b)

@deprecate csgintersection(a, b) a ∩ b
@deprecate csgintersection(a, b, tr) transform(a ∩ b, tr)

"""
    ∪(a::Union{CSGGenerator{T},ParametricSurface{T}}, b::Union{CSGGenerator{T},
        ParametricSurface{T}}) where {T<:Real}

Create a binary node in the CSG tree representing a union between `a` and `b`.

![Union Image](https://upload.wikimedia.org/wikipedia/commons/4/4a/Boolean_union.PNG)
"""
Base.:∪

function Base.:∪(a::CSGGenerator{T}, b::CSGGenerator{T}) where {T <: Real}
    return CSGGenerator{T}(
        (parenttransform) -> UnionNode(a(parenttransform), b(parenttransform)),
    )
end
Base.:∪(a::ParametricSurface, b::ParametricSurface) = leaf(a) ∪ leaf(b)
Base.:∪(a::ParametricSurface, b::CSGGenerator) = leaf(a) ∪ b
Base.:∪(a::CSGGenerator, b::ParametricSurface) = a ∪ leaf(b)

@deprecate csgunion(a, b) a ∪ b
@deprecate csgunion(a, b, tr) transform(a ∪ b, tr)

"""
    -(a::Union{CSGGenerator{T},ParametricSurface{T}}, b::Union{CSGGenerator{T},
        ParametricSurface{T}}) where {T<:Real}

Create a binary node in the CSG tree representing the difference of `a` and `b`,
essentially `a - b`.

![Difference Image](https://upload.wikimedia.org/wikipedia/commons/8/86/Boolean_difference.PNG)
"""
Base.:-

function Base.:-(a::CSGGenerator{T}, b::CSGGenerator{T}) where {T <: Real}
    return CSGGenerator{T}(
        (parenttransform) -> IntersectionNode(
            a(parenttransform),
            ComplementNode(b(parenttransform)),
        ),
    )
end
Base.:-(a::ParametricSurface, b::ParametricSurface) = leaf(a) - leaf(b)
Base.:-(a::ParametricSurface, b::CSGGenerator) = leaf(a) - b
Base.:-(a::CSGGenerator, b::ParametricSurface) = a - leaf(b)

@deprecate csgdifference(a, b) a - b
@deprecate csgdifference(a, b, tr) transform(a - b, tr)

"""
    evalcsg(
        a::Union{UnionNode{T},IntersectionNode{T},ComplementNode{T},LeafNode{T}},
        ray::AbstractRay{T,N},
        normalreverse::Bool = false
    )::Union{EmptyInterval{T},DisjointUnion{T},Interval{T}}

[TODO]
"""
function evalcsg end

function evalcsg(
    a::UnionNode{T},
    ray::AbstractRay{T, N},
    normalreverse::Bool = false,
) where {T <: Real, N}
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    l = evalcsg(a.leftchild, ray, normalreverse)
    if l isa Interval{T}
        if lower(l) isa RayOrigin{T} && upper(l) isa Infinity{T}
            return l
        end
    end
    r = evalcsg(a.rightchild, ray, normalreverse)
    if l isa EmptyInterval{T}
        return r
    end
    if r isa Interval{T}
        if lower(r) isa RayOrigin{T} && upper(r) isa Infinity{T}
            return r
        end
    end
    if r isa EmptyInterval{T}
        return l
    end
    return intervalunion(l, r)
end

function evalcsg(
    a::IntersectionNode{T},
    ray::AbstractRay{T, N},
    normalreverse::Bool = false,
) where {T <: Real, N}
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    l = evalcsg(a.leftchild, ray, normalreverse)
    if l isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        r = evalcsg(a.rightchild, ray, normalreverse)
        if r isa EmptyInterval{T}
            return EmptyInterval(T)
        else
            return intervalintersection(l, r)
        end
    end
end

function evalcsg(
    a::ComplementNode{T},
    ray::AbstractRay{T, N},
    normalreverse::Bool = false,
) where {T <: Real, N}
    return intervalcomplement(evalcsg(a.child, ray, !normalreverse))
end

function evalcsg(
    a::LeafNode{T},
    ray::AbstractRay{T, N},
    normalreverse::Bool = false,
) where {T <: Real, N}
    # the bounding box is in global coordinates so use the un-transformed ray
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    intscts = surfaceintersection(a.geometry, a.invtransform * ray)
    if !(intscts isa EmptyInterval{T})
        if normalreverse
            intscts = reversenormal(intscts)
        end
        if intscts isa Interval{T}
            return a.transform * intscts
        elseif intscts isa DisjointUnion{T}
            for i in eachindex(intervals(intscts))
                intervals(intscts)[i] = a.transform * intervals(intscts)[i]
            end
            return intscts
        else
            throw(
                ErrorException(
                    "EvalCSG error - returned type not an interval or disjoint union",
                ),
            )
        end
    else
        return EmptyInterval(T)
    end
end

"""
    surfaceintersection(obj::CSGTree{T}, r::AbstractRay{T,N})

Calculates the intersection of `r` with CSG object, `obj`.

Returns an [`EmptyInterval`](@ref) if there is no intersection, an
[`Interval`](@ref) if there is one or two interesections and a
[`DisjointUnion`](@ref) if there are more than two intersections.

The ray is intersected with the [`LeafNode`](@ref)s that make up the CSG object
and the resulting `Interval`s and `DisjointUnion`s are composed with the same
boolean operations to give a final result.
The ray is transformed by the inverse of the transform associated with the leaf
node to put it in _object space_ for
that node before the intersection is carried out, typically this _object space_
is centered at the origin, but may differ for each primitive.

Some intersections are culled without actually evaluating them by first checking
if the ray intersects the [`BoundingBox`](@ref) of each node in the
[`CSGTree`](@ref), this can substantially improve performance in some cases.
"""
function surfaceintersection(
    element::CSGTree{T},
    r::AbstractRay{T, N},
) where {T <: Real, N}
    return evalcsg(element, r)
end

point(::EmptyInterval) = nothing
normal(::EmptyInterval) = nothing
point(a::DisjointUnion) = point(closestintersection(a))
normal(a::DisjointUnion) = normal(closestintersection(a))
point(a::Interval) = point(closestintersection(a))
normal(a::Interval) = point(closestintersection(a))

"""
    onsurface(obj::CSGTree{T}, point::SVector{3,T}) -> Bool
    onsurface(obj::CSGTree{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _on_ the surface (i.e. shell) of `obj`.
"""
function onsurface(a::CSGTree{T}, x::T, y::T, z::T) where {T <: Real}
    return onsurface(a, SVector{3, T}(x, y, z))
end
function onsurface(a::ComplementNode{T}, point::SVector{3, T}) where {T <: Real}
    return onsurface(a.child, point)
end
function onsurface(
    a::IntersectionNode{T},
    point::SVector{3, T},
) where {T <: Real}
    return (
        (onsurface(a.leftchild, point) && inside(a.rightchild, point)) ||
        (onsurface(a.rightchild, point) && inside(a.leftchild, point))
    )
end
function onsurface(a::UnionNode{T}, point::SVector{3, T}) where {T <: Real}
    return (
        (onsurface(a.leftchild, point) && !inside(a.rightchild, point)) ||
        (onsurface(a.rightchild, point) && !inside(a.leftchild, point))
    )
end
function onsurface(a::LeafNode{T}, point::SVector{3, T}) where {T <: Real}
    return onsurface(a.geometry, a.invtransform * point)
end

"""
    inside(obj::CSGTree{T}, point::SVector{3,T}) -> Bool
    inside(obj::CSGTree{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _inside_ `obj`.
"""
function inside(a::CSGTree{T}, x::T, y::T, z::T) where {T <: Real}
    return inside(a, SVector{3, T}(x, y, z))
end
function inside(a::ComplementNode{T}, point::SVector{3, T}) where {T <: Real}
    return !inside(a.child, point)
end
function inside(a::IntersectionNode{T}, point::SVector{3, T}) where {T <: Real}
    return BoundingBoxes.inside(a.bbox, point) &&
           (inside(a.leftchild, point) && inside(a.rightchild, point))
end
function inside(a::UnionNode{T}, point::SVector{3, T}) where {T <: Real}
    return inside(a.bbox, point) &&
           (inside(a.leftchild, point) || inside(a.rightchild, point))
end
function inside(a::LeafNode{T}, point::SVector{3, T}) where {T <: Real}
    return BoundingBoxes.inside(a.bbox, point) &&
           inside(a.geometry, a.invtransform * point)
end

# ============================================================================ #

struct TrianglePool{T <: Real}
    allocated::Vector{Vector{Triangle{T}}}
    unallocated::Vector{Vector{Triangle{T}}}

    function TrianglePool{T}() where {T <: Real}
        return new{T}(
            Vector{Vector{Triangle{T}}}(),
            Vector{Vector{Triangle{T}}}(),
        )
    end
end

function allocate!(a::TrianglePool{T}) where {T <: Real}
    if length(a.unallocated) === 0
        #this will extend the array with a new Vector of Triangles.
        push!(a.unallocated, Vector{Triangle{T}}())
    end

    temp = pop!(a.unallocated)
    Base.empty!(temp) #resets count in array but doesn't reclaim space
    push!(a.allocated, temp)
    return temp
end

function empty!(a::TrianglePool{T}) where {T <: Real}
    for i in 1:length(a.allocated)
        push!(a.unallocated, pop!(a.allocated))
    end
end

# Allocate a zero length array which will be filled with Threads.nthreads()
# entries by the __init__ method for OpticSim module. Have to do this in the
# __init__ method because this captures the load time environment. const values
# are evaluated at precompile time and the number of threads in these two
# environments can be different.
const threadedtrianglepool = Vector{Dict{DataType, TrianglePool}}()

function newintrianglepool!(
    ::Type{T} = Float64,
    tid::Int = Threads.threadid(),
) where {T <: Real}
    if T ∉ keys(threadedtrianglepool[tid])
        # if the type of the interval pool has changed then we need to refill
        # it with the correct type
        threadedtrianglepool[tid][T] = TrianglePool{T}()
    end
    return allocate!(threadedtrianglepool[tid][T])
end

function emptytrianglepool!(
    ::Type{T} = Float64,
    tid::Int = Threads.threadid(),
) where {T}
    if T ∈ keys(threadedtrianglepool[tid])
        empty!(threadedtrianglepool[tid][T])
    end
end

# ============================================================================ #

# currrently we assume that no triangles span 2 sides of a shape, i.e. any
# triangle only crosses one shape boundary and so splitting is trivial
function uniontri!(
    csg::CSGTree{T},
    tri::Triangle{T},
    triangles::Vector{Triangle{T}},
    thisforcoplanar::Bool = false,
) where {T <: Real}
    v1 = vertex(tri, 1)
    v2 = vertex(tri, 2)
    v3 = vertex(tri, 3)
    in1 = inside(csg, v1)
    in2 = inside(csg, v2)
    in3 = inside(csg, v3)
    # second condition for if the vertex is in the plane of the other object,
    # in which case we want to keep one copy of the two coplanar surfaces
    if !thisforcoplanar
        # if we are NOT using this surface for coplanar faces then we want to
        # count the verts on the surface as being inside too in case we have
        # coplanar faces
        in1 |= onsurface(csg, v1)
        in2 |= onsurface(csg, v2)
        in3 |= onsurface(csg, v3)
    end
    if !(in1 || in2 || in3) || (
        thisforcoplanar &&
        onsurface(csg, v1) &&
        onsurface(csg, v2) &&
        onsurface(csg, v3)
    )
        # whole triangle is valid
        push!(triangles, tri)
    elseif !in1 && !in2 && in3
        splittri2out!(csg, v1, v2, v3, triangles)
    elseif !in2 && !in3 && in1
        splittri2out!(csg, v2, v3, v1, triangles)
    elseif !in1 && !in3 && in2
        splittri2out!(csg, v3, v1, v2, triangles)
    elseif !in1 && in2 && in3
        splittri1out!(csg, v1, v2, v3, triangles)
    elseif !in2 && in1 && in3
        splittri1out!(csg, v2, v3, v1, triangles)
    elseif !in3 && in1 && in2
        splittri1out!(csg, v3, v1, v2, triangles)
    end
end

function intersecttri!(
    csg::CSGTree{T},
    tri::Triangle{T},
    triangles::Vector{Triangle{T}},
    thisforcoplanar::Bool = false,
) where {T <: Real}
    v1 = vertex(tri, 1)
    v2 = vertex(tri, 2)
    v3 = vertex(tri, 3)
    in1 = inside(csg, v1)
    in2 = inside(csg, v2)
    in3 = inside(csg, v3)
    if thisforcoplanar
        # if we are allowing dupes then we want to count the verts on the
        # surface as being inside too in case we have coplanar faces
        in1 |= onsurface(csg, v1)
        in2 |= onsurface(csg, v2)
        in3 |= onsurface(csg, v3)
    end
    if (in1 && in2 && in3)
        # whole triangle is valid
        push!(triangles, tri)
    elseif in1 && in2 && !in3
        splittri2out!(csg, v1, v2, v3, triangles)
    elseif in2 && in3 && !in1
        splittri2out!(csg, v2, v3, v1, triangles)
    elseif in1 && in3 && !in2
        splittri2out!(csg, v3, v1, v2, triangles)
    elseif in1 && !in2 && !in3
        splittri1out!(csg, v1, v2, v3, triangles)
    elseif in2 && !in1 && !in3
        splittri1out!(csg, v2, v3, v1, triangles)
    elseif in3 && !in1 && !in2
        splittri1out!(csg, v3, v1, v2, triangles)
    end
end

function splittri1out!(
    csg::CSGTree{T},
    outv::SVector{3, T},
    inv1::SVector{3, T},
    inv2::SVector{3, T},
    triangles::Vector{Triangle{T}},
    recurse::Int = 0,
) where {T <: Real}
    if !validtri(outv, inv1, inv2) || recurse > VIS_RECURSION_LIMIT
        # if this triangle is invalid then just stop
        return
    end
    n1 = norm(outv - inv1)
    n2 = norm(outv - inv2)
    n3 = norm(inv1 - inv2)
    if min(n1, n2, n3) < MIN_VIS_TRI_SIZE
        # triangle is too small to bother splitting so just add it as is
        push!(triangles, Triangle(outv, inv1, inv2))
        return
    end
    n1 *= 1.0 + MESH_PRECISION
    n2 *= 1.0 + MESH_PRECISION
    int1 = closestintersection(evalcsg(csg, Ray(outv, inv1 - outv)), false)
    int2 = closestintersection(evalcsg(csg, Ray(outv, inv2 - outv)), false)
    if (
        !(int1 === nothing || int2 === nothing) &&
        validtri(outv, point(int1), point(int2)) &&
        α(int1) <= n1 &&
        α(int2) <= n2
    )
        mid = (inv1 + inv2) / 2
        test = closestintersection(evalcsg(csg, Ray(outv, mid - outv)), false)
        if test !== nothing
            if !isapprox(point(test), mid, atol = MESH_PRECISION)
                splittri1out!(
                    csg,
                    outv,
                    point(int1),
                    point(test),
                    triangles,
                    recurse + 1,
                )
                splittri1out!(
                    csg,
                    outv,
                    point(test),
                    point(int2),
                    triangles,
                    recurse + 1,
                )
            else
                push!(triangles, Triangle(outv, point(int1), point(int2)))
            end
        else
            push!(triangles, Triangle(outv, point(int1), point(int2)))
        end
    end
end

function splittri2out!(
    csg::CSGTree{T},
    outv1::SVector{3, T},
    outv2::SVector{3, T},
    inv::SVector{3, T},
    triangles::Vector{Triangle{T}},
    flip::Bool = false,
    recurse::Int = 0,
) where {T <: Real}
    if !validtri(outv1, outv2, inv) || recurse > VIS_RECURSION_LIMIT
        # if this triangle is invalid then just stop
        return
    end
    n1 = norm(outv1 - inv)
    n2 = norm(outv2 - inv)
    n3 = norm(outv1 - outv2)
    if min(n1, n2, n3) < MIN_VIS_TRI_SIZE
        # triangle is too small to bother splitting so just add it as is
        push!(triangles, Triangle(outv1, outv2, inv))
        return
    end
    n1 *= 1.0 + MESH_PRECISION
    n2 *= 1.0 + MESH_PRECISION
    int1 = closestintersection(evalcsg(csg, Ray(outv1, inv - outv1)), false)
    int2 = closestintersection(evalcsg(csg, Ray(outv2, inv - outv2)), false)
    if int2 !== nothing && α(int2) <= n2
        # int2 is valid
        if int1 !== nothing && α(int1) <= n1
            if isapprox(point(int1), point(int2), atol = MESH_PRECISION)
                # we just need one triangle in this case as inv lies on the
                # boundary of the shape
                if validtri(outv1, outv2, point(int2))
                    if flip
                        push!(triangles, Triangle(outv1, point(int2), outv2))
                    else
                        push!(triangles, Triangle(outv1, outv2, point(int2)))
                    end
                end
            else
                if flip
                    splittri1out!(
                        csg,
                        outv1,
                        point(int1),
                        point(int2),
                        triangles,
                    )
                else
                    splittri1out!(
                        csg,
                        outv1,
                        point(int2),
                        point(int1),
                        triangles,
                    )
                end
                splittri2out!(
                    csg,
                    outv2,
                    outv1,
                    point(int2),
                    triangles,
                    !flip,
                    recurse + 1,
                )
            end
        else
            if validtri(outv1, outv2, point(int2))
                if flip
                    push!(triangles, Triangle(outv1, point(int2), outv2))
                else
                    push!(triangles, Triangle(outv1, outv2, point(int2)))
                end
            end
        end
    elseif int1 !== nothing && α(int1) <= n1
        if validtri(outv1, outv2, point(int1))
            if flip
                push!(triangles, Triangle(outv1, point(int1), outv2))
            else
                push!(triangles, Triangle(outv1, outv2, point(int1)))
            end
        end
    end
end

# visualization functions from CSG objects, most actual work happens in the
# above functions
function makemeshi(
    a::UnionNode{T},
    subdivisions::Int = 30,
)::Vector{Triangle{T}} where {T <: Real}
    # do the two children in parallel
    th1 = Threads.@spawn makemesh(a.leftchild, subdivisions)
    th2 = Threads.@spawn makemesh(a.rightchild, subdivisions)
    ltris = (fetch(th1)::TriangleMesh{T}).triangles
    rtris = (fetch(th2)::TriangleMesh{T}).triangles
    triangles = newintrianglepool!(T)
    @inbounds for tri in ltris
        uniontri!(a.rightchild, tri, triangles, true)
    end
    @inbounds for tri in rtris
        uniontri!(a.leftchild, tri, triangles)
    end
    return triangles
end

function makemeshi(
    a::IntersectionNode{T},
    subdivisions::Int = 30,
)::Vector{Triangle{T}} where {T <: Real}
    # do the two children in parallel
    th1 = Threads.@spawn makemesh(a.leftchild, subdivisions)
    th2 = Threads.@spawn makemesh(a.rightchild, subdivisions)
    ltris = (fetch(th1)::TriangleMesh{T}).triangles
    rtris = (fetch(th2)::TriangleMesh{T}).triangles
    triangles = newintrianglepool!(T)
    @inbounds for tri in ltris
        intersecttri!(a.rightchild, tri, triangles, true)
    end
    @inbounds for tri in rtris
        intersecttri!(a.leftchild, tri, triangles)
    end
    return triangles
end

function makemeshi(
    a::ComplementNode{T},
    subdivisions::Int = 30,
)::Vector{Triangle{T}} where {T <: Real}
    child_triangles = makemeshi(a.child, subdivisions)
    # flip the normals
    triangles = newintrianglepool!(T)
    @inbounds for i in 1:length(child_triangles)
        tri = child_triangles[i]
        push!(
            triangles,
            Triangle(vertex(tri, 3), vertex(tri, 2), vertex(tri, 1)),
        )
    end
    return triangles
end

function makemeshi(
    a::LeafNode{T},
    subdivisions::Int = 30,
)::Vector{Triangle{T}} where {T <: Real}
    child_triangles = triangulate(a.geometry, subdivisions)
    triangles = newintrianglepool!(T)
    @inbounds for i in 1:length(child_triangles)
        push!(triangles, a.transform * child_triangles[i])
    end
    return triangles
end

function makemesh(
    c::CSGTree{T},
    subdivisions::Int = 30,
)::TriangleMesh{T} where {T <: Real}
    # need to copy with multithreading
    m = TriangleMesh(copy(makemeshi(c, subdivisions)))
    emptytrianglepool!(T)
    emptyintervalpool!(T)
    return m
end

"""
    triangulate(surf::ParametricSurface{S,N}, quads_per_row::Int, extensionu::Bool = false, extensionv::Bool = false, radialu::Bool = false, radialv::Bool = false)

Create an array of triangles representing the parametric surface where vertices are sampled on an even grid in UV space.
The surface can be extended by 1% in u and v separately, and specifying either u or v as being radial - i.e. detemining the radius on the surface e.g. rho for zernike - will result in that dimension being sampled using sqwrt so that area of triangles is uniform. The extension will also only apply to the maximum in this case.
"""
function triangulate(
    surf::ParametricSurface{T, N},
    subdivisons::Int,
    extensionu::Bool = false,
    extensionv::Bool = false,
    radialu::Bool = false,
    radialv::Bool = false,
) where {T, N}
    triangles = newintrianglepool!(T)
    (umin, umax), (vmin, vmax) = uvrange(surf)

    if extensionu
        if !radialu
            umin -= TRIANGULATION_EXTENSION
        end
        umax += TRIANGULATION_EXTENSION
    end

    if extensionv
        if !radialv
            vmin -= TRIANGULATION_EXTENSION
        end
        vmax += TRIANGULATION_EXTENSION
    end

    # if we are using this for intersection then we need to expand the triangles very slightly to
    # avoid misses due to floating point precision
    expansion = extensionu || extensionv ? TRIANGULATION_EXPANSION : zero(T)

    for ui in 0:(subdivisons - 1)
        for vj in 0:(subdivisons - 1)
            # this can be inefficient as we evaluate each point twice
            # in practice we usually use extension in which case no point is the same anyway
            u1 = T(ui / subdivisons)
            u2 = T((ui + 1) / subdivisons)
            v1 = T(vj / subdivisons)
            v2 = T((vj + 1) / subdivisons)

            if radialu
                u1 = sqrt(max(u1, zero(T)))
                u2 = sqrt(max(u2, zero(T)))
            end
            if radialv
                v1 = sqrt(max(v1, zero(T)))
                v2 = sqrt(max(v2, zero(T)))
            end

            u1 = (umax - umin) * u1 + umin - expansion
            u2 = (umax - umin) * u2 + umin + expansion
            v1 = (vmax - vmin) * v1 + vmin - expansion
            v2 = (vmax - vmin) * v2 + vmin + expansion

            p1 = point(surf, u1, v1)
            p2 = point(surf, u2, v1)
            p3 = point(surf, u2, v2)
            p4 = point(surf, u1, v2)

            if validtri(p1, p2, p3)
                push!(
                    triangles,
                    Triangle(
                        p1,
                        p2,
                        p3,
                        SVector{2, T}(u1, v1),
                        SVector{2, T}(u2, v1),
                        SVector{2, T}(u2, v2),
                    ),
                )
            end
            if validtri(p1, p3, p4)
                push!(
                    triangles,
                    Triangle(
                        p1,
                        p3,
                        p4,
                        SVector{2, T}(u1, v1),
                        SVector{2, T}(u2, v2),
                        SVector{2, T}(u1, v2),
                    ),
                )
            end
        end
    end

    return triangles
end

"""
    makemesh(object, subdivisions::Int = 30) -> TriangleMesh

Creates a [`TriangleMesh`](@ref) from an object, either a [`ParametricSurface`](@ref), [`CSGTree`](@ref) or certain surfaces (e.g. `Circle`, `Rectangle`).
This is used for visualization purposes only.
"""
function makemesh(
    surface::ParametricSurface{S, N},
    subdivisions::Int = 30,
)::TriangleMesh{S} where {S, N}
    m = TriangleMesh(triangulate(surface, subdivisions, false))
    emptytrianglepool!(S)
    return m
end

end # module


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using FactCheck
using Latte

function convolution_forward(filter::Array, bias::Array, input::Array,
                             output::Array, stride, pad, kernel)
    width, height, channel, num = size(input)
    width_out, height_out, n_filter, _ = size(output)

    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out,
        k = 1:channel, p = 1:kernel, q = 1:kernel

        in_y = (y-1) * stride - pad + p
        in_x = (x-1) * stride - pad + q
        if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
            filter_idx = ((k - 1) * kernel + (p - 1)) * kernel + q
            output[x, y, o, n] += input[in_x, in_y, k, n] *
                                  filter[filter_idx, o]
        end
    end

    # add bias
    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out
        output[x, y, o, n] += bias[o]
    end

    return output
end

function convolution_backward(filter::Array, bias::Array, input::Array,
                              top_diff::Array, stride, pad, kernel,
                              ∇filter::Array, ∇bias::Array, ∇input::Array)
    width, height, channels, num = size(input)
    width_out, height_out, n_filter, _ = size(top_diff)

    # ∇bias
    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out
        ∇bias[o] += top_diff[x, y, o, n]
    end

    # ∇filter and ∇input
    for n = 1:num, o = 1:n_filter, k = 1:channels, y = 1:height_out,
        x = 1:width_out, p = 1:kernel, q = 1:kernel

        in_y = (y-1) * stride - pad + p
        in_x = (x-1) * stride - pad + q
        if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
            filter_idx = ((k - 1) * kernel + (p - 1)) * kernel + q
            ∇filter[filter_idx, o] += top_diff[x,y,o,n] * input[in_x,in_y,k,n]
            ∇input[in_x,in_y,k,n] += top_diff[x,y,o,n] * filter[filter_idx,o]
        end
    end

    return (∇filter, ∇bias, ∇input)
end

pad = 1

srand(1234)
net = Net(8)
data,  data_value   = MemoryDataLayer(net, :data, (227, 227, 3))
label, label_value = MemoryDataLayer(net, :label, (1,))
data_value[:]  = rand(Float32, size(data_value)...) * 256
label_value[:] = map(floor, rand(Float32, size(label_value)...) * 10)
conv1        = ConvolutionLayer(:conv1, net, data, 10, 3, 1, pad)
conv2        = ConvolutionLayer(:conv2, net, conv1, 10, 3, 1, pad)
fc1          = InnerProductLayer(:fc1, net, conv2, 10)
loss         = SoftmaxLossLayer(:loss, net, fc1, label)

init(net)
ϵ = 1e-3

input    = get_buffer(net, :conv1value)
∇input   = get_buffer(net, :conv1∇)
filters  = get_buffer(net, :conv2weights)
∇filters = get_buffer(net, :conv2∇weights)
bias     = get_buffer(net, :conv2bias)
rand!(bias)
∇bias    = get_buffer(net, :conv2∇bias)
top_diff = get_buffer(net, :conv2∇)
# rand!(top_diff)
expected = zeros(get_buffer(net, :conv2value))
∇input_expected   = zeros(∇input)
∇filters_expected = zeros(size(∇filters)[1:end-1])
∇bias_expected    = zeros(size(∇bias)[1:end-1])

params = SolverParameters(
    lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
    mom_policy   = MomPolicy.Fixed(0.9),
    max_epoch    = 300,
    regu_coef    = .0005)
sgd = SGD(params)

facts("Testing Convolution Layer") do
    context("Forward") do
        forward(net; solver=sgd)
        convolution_forward(filters, bias, input, expected, 1, pad, 3)
        # @fact expected --> roughly(net.buffers[:conv2value])
        @fact all(-ϵ .< expected - get_buffer(net, :conv2value) .< ϵ) --> true
    end
    context("Backward") do
        clear_∇(net)
        backward(net)
        convolution_backward(filters, bias, input, top_diff, 1, pad, 3,
                             ∇filters_expected, ∇bias_expected,
                             ∇input_expected)
        @fact all(-ϵ .< ∇input - ∇input_expected .< ϵ) --> true
        @fact all(-ϵ .< sum(∇filters, 3)[:,:,1] - ∇filters_expected .< ϵ) --> true
        @fact all(-ϵ .< sum(∇bias, 3)[:,:,1] - ∇bias_expected .< ϵ) --> true
    end
end

FactCheck.exitstatus()

ext: jl


text:
using Test, Random
using Statistics, SparseArrays, LinearAlgebra
using Optim
using LogisticOptTools

const LOT = LogisticOptTools
const SVM_DATASET = joinpath(@__DIR__, "diabetes.txt")
Random.seed!(2713)

function fake_dataset(n_samples=100, n_features=10)
    X = randn(n_samples, n_features)
    w = randn(n_features)
    y = sign.(X * w)
    X .+= 0.8 * randn(n_samples, n_features) # add noise
    X .+= 100.0 # this makes it correlated by adding a constant term
    X = hcat(X, ones(n_samples, 1))
    return X, y
end

@testset "Utils functions" begin
    # Test expit
    @testset "Evaluation functions" begin
        for t in [-1.0, 0.0, 1.0]
            @test LOT.expit(t) == 1.0 / (1.0 + exp(-t))
            @test LOT.log1pexp(t) == -log(1.0 + exp(-t))
        end
    end
    @testset "Preprocessing" begin
        n, m = 100, 10
        dataset = randn(n, m)
        LOT.scale!(LOT.NormalScaler(), dataset)
        @test mean(dataset, dims=1) ≈ zeros(1, m) atol=1e-8
        @test std(dataset, dims=1) ≈ ones(1, m)

        y1 = [-2, 2, 2, -2]
        LOT.format_label!(y1)
        @test isequal(y1, [-1, 1, 1, -1])
        y2 = [-1, 3, 3]
        LOT.format_label!(y2)
        @test isequal(y2, [-1, 1, 1])
        y3 = [-1, 0, 0]
        LOT.format_label!(y3)
        @test isequal(y3, [-1, 1, 1])
    end
end

@testset "LIBSVM parser" begin
    # Load dataset from LIBSVM URL
    if !isfile(SVM_DATASET)
        download("$(LOT.LIBSVM_URL)/diabetes", SVM_DATASET)
    end
    # Parse diabetes either with Float32 or Float64
    for type_ in [Float32, Float64]
        dat = LOT.parse_libsvm(SVM_DATASET, type_)
        @test isa(dat, LOT.LibSVMData{Int, type_})
        X = LOT.to_dense(dat)
        @test isa(X, Array{type_, 2})
        Xs = LOT.to_sparse(dat)
        @test issparse(Xs)
    end
end

@testset "Primal model" begin
    # Test dataset
    svm_data = LOT.parse_libsvm(SVM_DATASET, Float64)
    X = LOT.to_dense(svm_data)
    y = svm_data.labels

    @testset "Constructor dense dataset" begin
        d1 = @inferred LOT.LogitData(X, y)
        @test sort(unique(d1.y)) == [-1.0, 1.0]
        d2 = @inferred LOT.LogitData(SVM_DATASET, scale_data=false)
        @test sort(unique(d2.y)) == [-1.0, 1.0]
        # Import should be the same when data are not scaled.
        @test hash(d1.X) == hash(d2.X)
        d3 = @inferred LOT.LogitData(SVM_DATASET, scale_data=true)
        @test sort(unique(d3.y)) == [-1.0, 1.0]
    end

    @testset "Constructor sparse dataset" begin
        d2 = @inferred LOT.SparseLogitData(SVM_DATASET, scale_data=false)
        @test sort(unique(d2.y)) == [-1.0, 1.0]
        d2 = @inferred LOT.SparseLogitData(SVM_DATASET, scale_data=true)
        @test sort(unique(d2.y)) == [-1.0, 1.0]
    end

    @testset "Primal logistic regression" begin
        for Dataset in [LOT.LogitData, LOT.SparseLogitData]
            data = Dataset(SVM_DATASET)
            # Input dimension
            p = LOT.nfeatures(data)
            @test isa(p, Int)
            @test p == size(data.X, 2)
            # Initial point
            x0 = zeros(p)
            # vector for hess-vec product
            vec = zeros(p)
            # Allocate vectors for gradient, Hessian and diag-Hessian
            g = zeros(p)
            hess = zeros(div(p * (p+1), 2))
            diagh = zeros(p)
            @test LOT.loss(x0, data) ≈ -LOT.log1pexp(0.0)
            LOT.gradient!(g, x0, data)
            LOT.hessvec!(g, x0, vec, data)
            LOT.diaghess!(diagh, x0, data)
            if isa(data, LOT.LogitData)
                LOT.hessian!(hess, x0, data)
            end
        end
    end

    @testset "Intercept" begin
        data = LOT.LogitData(SVM_DATASET)
        # +1 accounts for intercept
        p = LOT.nfeatures(data) + 1
        # Initial point
        x0 = zeros(p)
        # vector for hess-vec product
        vec = zeros(p)
        # Allocate vectors for gradient, Hessian and diag-Hessian
        g = zeros(p)
        hess = zeros(div(p * (p+1), 2))
        diagh = zeros(p)
        fit_intercept = true
        LOT.gradient!(g, x0, data, fit_intercept)
        LOT.hessvec!(g, x0, vec, data, fit_intercept)
        LOT.diaghess!(diagh, x0, data, fit_intercept)
        LOT.hessian!(hess, x0, data, fit_intercept)
    end

    @testset "Penalty" begin
        # Penalty parameter
        p = size(X, 2)
        x0 = zeros(p)
        g = zeros(p)
        hess = zeros(div(p * (p+1), 2))
        diagh = zeros(p)
        fill!(g, 0.0)
        fill!(hess, 0.0)
        fill!(diagh, 0.0)
        λ = 1.0
        for penalty in [LOT.L2Penalty(λ),
                        LOT.L1Penalty(λ),
                        LOT.LinearizedL1Penalty(λ),
                        LOT.L0Penalty(1.0)]
            @test penalty(x0) == 0.0
            LOT.gradient!(g, x0, penalty)
            @test all(g .== 0.0)
            LOT.hessian!(hess, x0, penalty)
            LOT.hessvec!(g, x0, vec, penalty)
            LOT.diaghess!(diagh, x0, penalty)
            @test all(diagh .== 2.0)
        end
    end

    @testset "Fitting logistic regression (primal)" begin
        # Load dataset
        svm_data = LOT.parse_libsvm(SVM_DATASET, Float64)
        X1 = LOT.to_dense(svm_data)
        X2 = LOT.to_sparse(svm_data)
        y = svm_data.labels
        # Solution when data are scaled
        solution = 0.6084979240046
        options = Optim.Options(iterations=250, g_tol=1e-5)
        algo = BFGS()
        # No intercept
        for (X, glm) in zip([X1, X2], [LOT.LogitData, LOT.SparseLogitData])
            data = glm(X, y)
            @test LOT.ndata(data) == length(y)
            p = LOT.nfeatures(data)
            @test p == size(X, 2)
            # Build callbacks
            logreg = LOT.LogisticRegressor(X, y)
            f, g!, _, _ = LOT.generate_callbacks(logreg)
            x0 = zeros(p)
            res_joptim = Optim.optimize(f, g!, x0, algo, options)
            @test res_joptim.minimum ≈ solution
            @test Optim.converged(res_joptim)
            # Test dedicated shortcut function
            res_joptim = Optim.optimize(logreg, x0, algo, options)
            @test res_joptim.minimum ≈ solution
            @test Optim.converged(res_joptim)
        end

        # Intercept
        solution_intercept = 0.47099308415704666
        for (X, glm) in zip([X1, X2], [LOT.LogitData, LOT.SparseLogitData])
            data = glm(X, y)
            p = LOT.nfeatures(data)
            # Build callbacks
            logreg = LOT.LogisticRegressor(X, y, fit_intercept=true)
            f, g!, _, _ = LOT.generate_callbacks(logreg)
            x0 = zeros(p + 1)
            res_joptim = Optim.optimize(f, g!, x0, algo, options)
            @test res_joptim.minimum ≈ solution_intercept
            @test Optim.converged(res_joptim)
            # Test dedicated shortcut function
            res_joptim = Optim.optimize(logreg, x0, algo, options)
            @test res_joptim.minimum ≈ solution_intercept
            @test Optim.converged(res_joptim)
        end

        @testset "LogisticOptimizer utilities" begin
            # Test fit! function
            for (bias, sol) in zip([false, true], [solution, solution_intercept])
                logopt = LOT.LogisticOptimizer(algo=algo, options=options,
                                               fit_intercept=bias)
                p = size(X, 2) + bias
                res = LOT.fit!(logopt, X1, y, zeros(p))
                @test Optim.converged(res)
                @test res.minimum ≈ sol
            end
        end
    end
end

@testset "Test optimization of L2 penalty" begin
    dataset = LOT.LogitData(SVM_DATASET, scale_data=true)
    kfold = 5
    shuffle = false
    inner_algo = LBFGS()
    outer_algo = BFGS()
    options = Optim.Options()
    penopt = LOT.L2PenaltyOptimizer(kfold, shuffle, false,
                                    outer_algo, inner_algo, options)
    l♯, γ♯ = LOT.fit!(penopt, dataset.X, dataset.y, 1.0)
    @test l♯ ≈ 0.5536870145713175 atol=1e-4
    @test γ♯ ≈ 0.015463419098934016 atol=1e-4
end

@testset "Test dual model" begin
    svm_data = LOT.parse_libsvm(SVM_DATASET, Float64)
    X = LOT.to_dense(svm_data)
    y = svm_data.labels
    data = LOT.DualLogitData(X, y)
    n = LOT.ndata(data)
    penalty = LOT.L2Penalty(0.0)
    @testset "Dual Logit" begin
        λ = zeros(n)
        @test LOT.loss(λ, data) == 0.0
        λ = -data.y
        @test LOT.loss(λ, data) == 0.0
        g = zeros(n)
        LOT.gradient!(g, λ, data)
    end
    @testset "Dual GLM" begin
        λ = zeros(n)
        g = zeros(n)
        glm = LOT.DualLogisticRegressor(data, penalty)
        @test LOT.loss(λ, glm) == 0.0
        LOT.gradient!(g, λ, glm)

        glm = LOT.DualLogisticRegressor(X, y)
        @test LOT.loss(λ, glm) == 0.0
        LOT.gradient!(g, λ, glm)
    end
    @testset "Fitting (dual)" begin
        glm = LOT.DualLogisticRegressor(data, penalty)
        f, gradient!, _, _ = LOT.generate_callbacks(glm)
        algo = BFGS()
        options = Optim.Options(iterations=250, g_tol=1e-5)
        lower = LOT.lowerbound(data)
        upper = LOT.upperbound(data)
        x0 = 0.5 * (lower .+ upper)
        res_joptim = Optim.optimize(f, gradient!, lower, upper,
                                    x0, Fminbox(algo), options)
        @test Optim.converged(res_joptim)
    end
end


text:
const BaseRules = OrderedRuleSet(
    "Base",
    ("LowCounts", r -> get(r, :COUNTS, 2001.0) < 2000.0),
    (
        "Maraging",
        r -> (r.FIRSTELM == n"Fe") &&
                 (r.NI > 5) && (r.CO > 2) && sum(r, :FE, :NI, :CO, :MO) > 80.0,
    ),
    ("Quartz-like", r -> r.SI > 80),
    (
        "Biotite",
        r -> (r.K > 5) &&
                 (r.MG + r.FE > 20) &&
                 (r.AL > 10) &&
                 (r.SI > 10) && (r.K + r.MG + r.FE + r.AL + r.SI > 80) && (r.O > 5),
    ),
    (
        "Plagioclaise",
        r -> (r.CA + r.NA > 20) &&
                 (r.AL > 10) && (r.SI > 10) && (r.NA + r.CA + r.AL + r.SI > 80) && (r.O > 5),
    ),
    (
        "Othoclaise",
        r -> (r.K > 5) &&
                 (r.AL > 20) && (r.SI > 20) && (r.O > 10) && (r.K + r.AL + r.SI > 80),
    ),
    ("Dolomite", r -> (r.CA > 40) && (r.MG > 10) && (r.CA + r.MG > 70)),
    ("Olivine", r -> (r.MG + r.FE > 40) && (r.SI > 20) && (r.MG + r.FE + r.SI > 80)),
    (
        "Feldspar",
        r -> (r.K + r.NA + r.CA > 20) &&
                 (r.AL > 10) && (r.SI > 20) && (r.K + r.NA + r.CA + r.AL + r.SI > 80),
    ),
    ("Willemite", r -> (r.ZN > 20) && (r.SI > 20) && (r.ZN + r.SI > 80)),
    (
        "Garnet",
        r -> (r.MG + r.FE + r.MN + r.CA > 10) &&
                 (r.AL > 5) && (r.SI > 10) && (r.MG + r.FE + r.MN + r.CA + r.AL + r.SI > 80),
    ),
    (
        "Al-Cu alloy (2XXX)",
        r -> (r.AL > 60) && (r.SECONDELM == n"Cu") && (r.CU > 3) && (r.AL + r.CU > 90),
    ),
    (
        "Al-Mn alloy (3XXX)",
        r -> (r.AL > 60) && (r.SECONDELM == n"Mn") && (r.MN > 3) && (r.AL + r.MN > 90),
    ),
    (
        "Al-Si alloy (4XXX)",
        r -> (r.AL > 80) && (r.SECONDELM == n"Si") && (r.SI > 1) && (r.AL + r.SI > 90),
    ),
    (
        "Al-Mg alloy (5XXX)",
        r -> (r.FIRSTELM == n"Al") && (r.SECONDELM == n"Mg") && (r.MG > 5) && (r.AL + r.MG > 80),
    ),
    (
        "Al-Si-Mg (6XXX)",
        r -> (r.AL > 90) &&
                 (r.SI > 0) && (r.MG > 0) && (r.AL + r.SI + r.MG + r.MN + r.CU + r.FE > 98),
    ),
    (
        "Al-Zn alloy (7XXX)",
        r -> (r.AL > 80) && (r.SECONDELM == n"Zn") && (r.ZN > 2) && (r.AL + r.ZN > 90),
    ),
    ("Al-Fe alloy", r -> (r.AL > 80) && (r.SECONDELM == n"Fe") && (r.FE > 2) && (r.AL + r.FE > 90)),
    ("Aluminosilicate", r -> (r.AL > 30) && (r.SI > 10) && (r.AL + r.SI > 80)),
    ("Zircon", r -> (r.ZR > 20) && (r.ZR + r.SI > 80)),
    (
        "Elgiloy",
        r -> (r.CO > 20) &&
                 (r.CR > 10) &&
                 (r.NI > 5) &&
                 (r.FE > 5) &&
                 ((r.MO > 2) || (r.S > 2)) &&
                 (r.CO + r.CR + r.NI + r.FE + r.MO + r.MN + r.S > 80),
    ),
    ("Aluminum", r -> r.AL > 80),
    ("Fe-Al alloy", r -> (r.AL > 20) && (r.FE > 20) && (r.AL + r.FE > 80)),
    ("Calcium silicate", r -> (r.CA > 40) && (r.SI > 10) && (r.CA + r.SI > 80)),
    ("Fe-Cr stainless", r -> (r.FE > 50) && (r.CR > 5) && (r.NI < 3) && (r.FE + r.CR > 90)),
    (
        "Fe-Cr-Ni Stainless",
        r -> (r.FIRSTELM == n"Fe") && (r.CR > 5) && (r.NI > 5) && (r.FE + r.CR + r.NI > 80),
    ),
    ("Fe-Ni stainless", r -> (r.FE > 50) && (r.NI > 10) && (r.FE + r.NI > 90)),
    ("Iron oxide", r -> (r.FE > 80) && (r.O > 10)),
    ("Iron", r -> r.FE > 80),
    ("Calcite", r -> (r.CA > 80) && (r.O > 5)),
    ("Zinc", r -> r.ZN > 80),
    ("Silicate", r -> r.FIRSTELM == n"Si"),
    (
        "Anglesite",
        r -> (r.FIRSTELM == n"Pb") && (r.SECONDELM == n"S") && (r.PB + r.S > 80) && (r.O > 5),
    ),
    ("Hashemite", r -> (r.BA > 20) && (r.CR > 20) && (r.BA + r.CR > 80) && (r.O > 5)),
    ("Barite", r -> (r.BA > 40) && (r.S > 10) && (r.BA + r.S > 80) && (r.O > 5)),
    ("Gypsum", r -> (r.CA > 20) && (r.S > 20) && (r.CA + r.S > 80) && (r.O > 5)),
    ("Lead", r -> (r.PB > 60) && (r.PB + r.S > 80)),
    ("Salt", r -> (r.K + r.NA > 20) && (r.CL > 20) && (r.K + r.NA + r.CL > 80)),
    ("Bismuth", r -> (r.BI > 80)),
    ("Brass", r -> (r.CU > 20) && (r.ZN > 20) && (r.CU + r.ZN > 80)),
    ("Bronze", r -> (r.CU > 20) && (r.SN > 20) && (r.CU + r.SN > 80)),
    ("Cupronickel", r -> (r.CU > 30) && (r.NI > 10) && (r.CU + r.NI > 80)),
    ("Sodium sulfide", r -> (r.NA > 20) && (r.S > 10) && (r.NA + r.S > 80)),
    ("Wurtzite (ZnS)", r -> (r.ZN > 20) && (r.S > 10) && (r.ZN + r.S > 80)),
    ("Pyrite (FeS2)", r -> (r.FE > 20) && (r.S > 20) && (r.FE + r.S > 80)),
    ("Silver sulfide", r -> (r.AG > 20) && (r.S > 10) && (r.AG + r.S > 80)),
    (
        "Sn62 solder",
        r -> (r.SN > 30) && (r.PB > 10) && (r.AG > 1) && (r.SN + r.PB + r.AG > 80),
    ),
    (
        "SACZ solder",
        r -> (r.SN + r.AG + r.CU + r.ZN > 80) &&
                 (r.SN > 10) && (r.CU > 10) && (r.AG > 10) && (r.ZN > 2),
    ),
    ("Lead solder", r -> (r.PB > 20) && (r.SN > 20) && (r.PB + r.SN > 80)),
    (
        "SAC solder",
        r -> (r.SN + r.AG + r.CU > 80) && (r.SN > 10) && (r.CU > 10) && (r.AG > 10),
    ),
    ("Titanium", r -> (r.TI > 80)),
    ("Calcium phosphate", r -> (r.CA > 40) && (r.P > 10) && (r.CA + r.P > 80)),
    ("Calcium fluoride", r -> (r.CA > 40) && (r.F > 10) && (r.CA + r.F > 80)),
    ("Iron-bearing", r -> (r.FE > 50)),
    ("Ca-bearing", r -> r.CA > 50),
    ("Other Silicate", r -> (r.SI > 20) && (r.AL + r.FE + r.CA + r.SI > 80)),
    ("Gold-Silver alloy", r -> (r.FIRSTELM == n"Au") && (r.SECONDELM == n"Ag") && (r.AU + r.AG > 80)),
    ("Gold - other", r -> r.AU > 80),
    ("Tin", r -> r.SN > 80),
    ("Celestine", r -> (r.SR > 40) && (r.S > 10) && (r.SR + r.S > 80) && (r.O > 5)),
    ("Barite-Celestine", r -> (r.BA + r.SR > 40) && (r.BA + r.SR + r.S > 80) && (r.O > 5)),
    ("Cu-rich", r -> r.CU > 80),
    ("Cr-bearing", r -> r.CR > 50),
    ("Silver sulfide", r -> (r.AG > 40) && (r.S > 10) && (r.AG + r.S > 80)),
    ("Silver", r -> r.AG > 80),
    ("MoS2", r -> r.MO > 40 && r.S > 30 && (r.MO + r.S > 90)),
    ("Monazite", r -> sum(r, :LA, :CE, :ND, :TH, :Y) > 30 && (r.P > 10) && sum(r, :LA, :CE, :ND, :TH, :Y, :P, :SI) > 80 ),
    ("Lanthanide", r -> sum(r, :LA, :CE, :ND, :Y) > 60),
    ("Silver chloride", r -> (r.AG > 40) && (r.CL > 10)),
    ("Magnesium oxide", r -> (r.MG > 80) && (r.O > 10)),
    ("Nickel", r -> r.NI > 80),
    ("Ca-Fe silicate", r -> (r.CA > 20) && (r.FE > 20) && (r.SI > 5)),
    ("Ca-Ti silicate", r -> (r.CA > 20) && (r.TI > 20) && (r.SI > 10)),
    ("Al-Zr-Cl", r -> (r.AL > 20) && (r.ZR > 20) && (r.CL > 5)),
    ("Chlorine-rich", r -> r.CL > 80),
    ("Au-Cu alloy", r -> (r.AU > 40) && (r.CU > 10) && (r.AU + r.CU + r.NI > 80)),
    ("Uranium-bearing", r -> r.U > 20),
    ("Thorium-bearing", r -> r.TH > 20),
    ("Platinum", r -> r.PT > 80)
)

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
abstract type PlotObject end

abstract type Axes <: PlotObject end

abstract type Plot <: PlotObject end

const FIGUREARGS = (:width, :height, :subplot, :currentaxes)

struct Figure
    axes::Dict{UInt8,Axes}
    args::Dict{Symbol,Any}
    function Figure(; kwargs...)
        fig = new(Dict{UInt8,Axes}(), Dict{Symbol,Any}())
        for (key, value) in kwargs
            key ∈ FIGUREARGS && (fig[key] = value)
        end
        get!(fig.args, :width, 604.72440945)
        get!(fig.args, :height, 427.60473067)
        get!(fig.args, :currentaxes, 1)
        get!(fig.args, :subplot, (1, 1))
        fig
    end
end

const CURRENTFIGURE = Ref(Figure())

Base.showable(::MIME"image/svg+xml", _::SVGPlots.PlotObject) = true

function Base.show(io::IO, mime::MIME"image/svg+xml", _::SVGPlots.PlotObject)
    show(io, mime, gcf())
end

Base.showable(::MIME"image/svg+xml", _::SVGPlots.Figure) = true

function Base.show(io::IO, mime::MIME"image/svg+xml", fig::SVGPlots.Figure)
    show(io, mime, Drawing(fig))
end

function Base.getindex(object::Union{PlotObject,Figure}, key::Symbol)
    return object.args[key]
end

function Base.setindex!(object::Union{PlotObject,Figure}, v, key::Symbol)
    object.args[key] = v
end

function savefig(filename::AbstractString, fig::Figure = gcf())
    write(filename, Drawing(fig))
end

function positionaxes(fig::Figure = gcf())
    xdim, ydim = fig[:subplot]
    ca = fig[:currentaxes]
    w = 1 / xdim
    h = 1 / ydim
    x, y = divrem(ca - 1, ydim)
    x * w, y * h, w, h
end

const CARTESIANARGS = (:hold, :xlim, :ylim, :xscale, :yscale, :position, :title, :xlabel, :ylabel, :legend)

struct Cartesian <: Axes
    plots::Vector{Plot}
    args::Dict{Symbol,Any}
    function Cartesian(; kwargs...)
        axes = new(Plot[], Dict{Symbol,Any}())
        for (key, value) in kwargs
            key ∈ CARTESIANARGS && (axes[key] = value)
        end
        get!(axes.args, :hold, :false)
        get!(axes.args, :xscale, :linear)
        get!(axes.args, :yscale, :linear)
        get!(axes.args, :position, positionaxes())
        axes
    end
end

const LINEARGS = (:linetype,)
const LINETYPES = Dict(
    :solid => [0],
    :dotted => [1 2],
    :dashed => [4 2],
    :dashdotted => [4 2 1 2],
)

struct Line <: Plot
    data::Array{Float64,2}
    args::Dict{Symbol,Any}
    function Line(x::AbstractVector, y::AbstractVector; kwargs...)
        length(x) != length(y) && error("Length not equal!")
        data = [float(x[:]) float(y[:])]
        line = new(data, Dict{Symbol,Any}())
        for (key, value) in kwargs
            key ∈ LINEARGS && (line[key] = value)
        end
        get!(line.args, :linetype, :solid)
        line
    end
end

gcf() = CURRENTFIGURE[]

function figure(fig::Union{Figure,Nothing} = nothing; kwargs...)
    fig == nothing && (fig = Figure(; kwargs...))
    CURRENTFIGURE[] = fig
end

function gca(fig::Figure = gcf())
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    fig.axes[ca]
end

function hold(val::Bool, fig::Figure = gcf(); kwargs...)
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    axes[:hold] = val
    axes
end

function subplot(xdim::Int, ydim::Int, ca::Int, fig::Figure = gcf())
    (xdim, ydim) != fig[:subplot] && empty!(fig.axes)
    fig[:currentaxes] = ca
    fig[:subplot] = (xdim, ydim)
    fig
end

function title(txt::String, fig::Figure = gcf(); kwargs...)
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    axes[:title] = txt
    axes
end

function xlabel(txt::String, fig::Figure = gcf(); kwargs...)
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    axes[:xlabel] = txt
    axes
end

function ylabel(txt::String, fig::Figure = gcf(); kwargs...)
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    axes[:ylabel] = txt
    axes
end

function legend(legs::Vector{S}, fig::Figure = gcf(); kwargs...) where S<:AbstractString
    ca = fig[:currentaxes]
    ca ∉ keys(fig.axes) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    axes[:legend] = legs
    axes
end

function plot(x::AbstractVector, y::AbstractVector, fig::Figure = gcf(); kwargs...)
    ca = fig[:currentaxes]
    ! (ca ∈ keys(fig.axes) && fig.axes[ca][:hold]) && (fig.axes[ca] = Cartesian(; kwargs...))
    axes = fig.axes[ca]
    line = Line(x, y; kwargs...)
    push!(axes.plots, line)
    line
end

const COLORS = (
    "rgb(0%, 44.7%, 74.1%)",
    "rgb(85%, 32.5%, 9.8%)",
    "rgb(92.9%, 69.4%, 12.5%)",
    "rgb(49.4%, 18.4%, 55.6%)",
    "rgb(46.6%, 67.4%, 18.8%)",
    "rgb(30.1%, 74.5%, 93.3%)",
    "rgb(63.5%, 7.8%, 18.4%)"
)

function Drawing(fig::Figure)
    width = fig[:width]
    height = fig[:height]
    NativeSVG.Drawing(
        width = width,
        height = height,
        xmlns!xlink = "http://www.w3.org/1999/xlink",
    ) do
        defs() do
            clipPath(id = "clipwindow") do
                rect(x = 0, y = 0, width = "100%", height = "100%")
            end
        end
        rect(x = 0, y = 0, width = "100%", height = "100%", fill = "white")
        g(clip_path = "url(#clipwindow)") do
            for axes in values(fig.axes)
                draw(axes, fig)
            end
        end
    end
end

function adjustlimits(amin::Number, amax::Number)
    a = log10(amax - amin)
    quotient, remainder = divrem(a, 1)
    remainder *= remainder
    exponent = floor(quotient)
    quotient - exponent > 0.5 && (exponent += 1)
    remainder < 0.5 && (exponent -= 1)
    scale = 10^(-exponent)
    floor(scale * amin) / scale, ceil(scale * amax) / scale
end

function tick(amin::Number, amax::Number)
    exponent = log10(amax - amin)
    intpart = div(exponent, 1)
    factor = 10^(exponent - intpart)
    tickunit = if factor > 5; 2
        elseif factor > 2.5; 1
        elseif factor > 1; 0.5
        elseif factor > 0.5; 0.2
        elseif factor > 0.25; 0.1
        else; 0.05
    end
    tickunit * 10^intpart
end

function draw(axes::Cartesian, fig::Figure)
    width = fig[:width]
    height = fig[:height]
    pos = axes[:position]
    x = pos[1] * width
    y = pos[2] * height
    w = pos[3] * width
    h = pos[4] * height
    fontsize = 12
    :title ∈ keys(axes.args) && text(
        x = x + 0.125 * w + 0.4 * w,
        y = y + 0.04 * h,
        text_anchor = :middle,
        dominant_baseline = :middle,
        font_size = fontsize,
        vector_effect="non-scaling-stroke"
    ) do
        str(axes[:title])
    end
    :xlabel ∈ keys(axes.args) && text(
        x = x + 0.125 * w + 0.4 * w,
        y = y + h - 0.02 * h,
        text_anchor = :middle,
        dominant_baseline = :middle,
        font_size = fontsize,
        font_style = :italic,
        vector_effect="non-scaling-stroke"
    ) do
        str(axes[:xlabel])
    end
    :ylabel ∈ keys(axes.args) && text(
        x = x + 0.02 * w,
        y = y + 0.075 * h + 0.4 * h,
        text_anchor = :middle,
        dominant_baseline = :middle,
        font_size = fontsize,
        font_style = :italic,
        rotate = -90,
        vector_effect="non-scaling-stroke"
    ) do
        str(axes[:ylabel])
    end
    (xmin, xmax) = if :xlim ∈ keys(axes.args)
        axes[:xlim]
    else
        xmin, xmax = 0, 0
        for plt in axes.plots
            xlim = extrema(plt.data[:, 1])
            xlim[1] < xmin && (xmin = xlim[1])
            xlim[2] > xmax && (xmax = xlim[2])
        end
        xmin, xmax = adjustlimits(xmin, xmax)
    end
    (ymin, ymax) = if :ylim ∈ keys(axes.args)
        axes[:ylim]
    else
        ymin, ymax = 0, 0
        for plt in axes.plots
            ylim = extrema(plt.data[:, 2])
            ylim[1] < ymin && (ymin = ylim[1])
            ylim[2] > ymax && (ymax = ylim[2])
        end
        ymin, ymax = adjustlimits(ymin, ymax)
    end
    a = 0.8 * w / (xmax - xmin)
    b = 0.125 * w + x - a * xmin
    c = 0.8 * h / (ymin - ymax)
    d = 0.075 * h + y - c * ymax
    xtick = tick(xmin, xmax) / 5
    tmp = xmin / xtick
    intpart = div(tmp, 1)
    i = intpart == tmp ? tmp : floor(intpart) + 1
    xi = i * xtick
    while xi <= xmax
        color, len = if rem(i, 5) == 0
            text(x = a * xi + b,
                y = c * ymin + d + 0.02 * h,
                text_anchor = :middle,
                dominant_baseline = :middle,
                font_size = fontsize,
                vector_effect="non-scaling-stroke"
            ) do
                str(string(xi == div(xi, 1) ? Int(xi) : xi))
            end
            :lightgrey, 0.02 * h
        else
            :whitesmoke, 0.01 * h
        end
        line(x1 = a * xi + b, y1 = c * ymin + d, x2 = a * xi + b, y2 = c * ymax + d, stroke = color)
        line(x1 = a * xi + b, y1 = c * ymin + d, x2 = a * xi + b, y2 = c * ymin + d - len, stroke = :black)
        line(x1 = a * xi + b, y1 = c * ymax + d, x2 = a * xi + b, y2 = c * ymax + d + len, stroke = :black)
        i += 1
        xi = i * xtick
    end
    ytick = tick(ymin, ymax) / 5
    tmp = ymin / ytick
    intpart = div(tmp, 1)
    i = intpart == tmp ? tmp : floor(intpart) + 1
    yi = i * ytick
    while yi <= ymax
        color, len = if rem(i, 5) == 0
            text(x = a * xmin + b - 0.02 * w,
                y = c * yi + d,
                text_anchor = :end,
                dominant_baseline = :middle,
                font_size = fontsize,
                vector_effect="non-scaling-stroke"
            ) do
                str(string(yi == div(yi, 1) ? Int(yi) : yi))
            end
            :lightgrey, 0.02 * h
        else
            :whitesmoke, 0.01 * h
        end
        line(x1 = a * xmin + b, y1 = c * yi + d, x2 = a * xmax + b, y2 = c * yi + d, stroke = color, vector_effect="non-scaling-stroke")
        line(x1 = a * xmin + b, y1 = c * yi + d, x2 = a * xmin + b + len, y2 = c * yi + d, stroke = :black, vector_effect="non-scaling-stroke")
        line(x1 = a * xmax + b, y1 = c * yi + d, x2 = a * xmax + b - len, y2 = c * yi + d, stroke = :black, vector_effect="non-scaling-stroke")
        i += 1
        yi = i * ytick
    end
        defs() do
            clipPath(id = "clipplot") do
                rect(x = 0, y = 0, width = "100%", height = "100%")
            end
        end
        for (i, plt) in enumerate(axes.plots)
            xx = a .* plt.data[:, 1] .+ b
            yy = c .* plt.data[:, 2] .+ d
            plt[:linetype] != :none &&
                polyline(
                    clipPath = "url(#clipplot)",
                    fill = "none",
                    stroke = COLORS[(i - 1) % 7 + 1],
                    stroke_dasharray = mapreduce(
                        e -> "$e ",
                        *,
                        LINETYPES[plt[:linetype]]
                    ),
                    points = mapreduce(
                        p -> "$(p[1]),$(p[2]) ",
                        *,
                        zip(xx, yy)
                    ),
                    vector_effect="non-scaling-stroke"
                )
        end
    rect(
        x = x + 0.125 * w,
        y = y + 0.075 * h,
        width = 0.8 * w,
        height = 0.8 * h,
        fill = "none",
        stroke = "black",
        vector_effect="non-scaling-stroke"
    )
end

max_stars_repo_name: JuliaTagBot/SVGPlots

max_issues_repo_name: JuliaTagBot/SVGPlots

max_forks_repo_name: JuliaTagBot/SVGPlots


text:
include("header.jl")
using Statistics, LinearAlgebra

sizes = [((2,4,3),(4,1,3)),((2,4,5),(4,8,5)),((2,8,4,3),(8,2,4,3))]
@testset "bmm" begin
    for t in (Float32, Float64)
        for s in sizes
            a = t(0.1)*randn(t, s[1]...)
            b = t(0.1)*randn(t, s[2]...)
            bmmul1(w)=bmm(w[1],w[2])
            bmmul2(w)=bmm(w[1],w[2]; transA=true)
            bmmul3(w)=bmm(w[1],w[2]; transB=true)
            bmmul4(w)=bmm(w[1],w[2]; transA=true, transB=true)
            pm(w) = ndims(w)==3 ? permutedims(w, (2,1,3)) : permutedims(w, (2,1,3,4))
            w = [a,b]
            @test gradcheck(bmmul1, w)
            w = [pm(a),b]
            @test gradcheck(bmmul2, w)
            w = [a,pm(b)]
            @test gradcheck(bmmul3, w)
            w = [pm(a),pm(b)]
            @test gradcheck(bmmul4, w)
            if gpu() >= 0
                c = bmm(a, b)
                ka = KnetArray(a)
                kb = KnetArray(b)
                kc = bmm(ka, kb)
                @test isapprox(c, Array(kc))
                w = [ka,kb]
                @test gradcheck(bmmul1, w)
            end
        end
    end
    if gpu() >= 0
        # Issue #495: transpose support
        ϵ =  1e-9
        A  = KnetArray(rand(3,4,10))
        At = permutedims(A,(2,1,3))
        B  = KnetArray(rand(4,5,10))
        Bt = permutedims(B,(2,1,3))
        C =  bmm(A,B)
        @test bmm(A,Bt;transB=true) ≈ C
        @test bmm(At,B;transA=true) ≈ C
        @test bmm(At,Bt;transA=true,transB=true) ≈ C
    end
end #end of testset

# suppress the return
nothing


text:
module TestBatchLenses

include("preamble.jl")
using Kaleido: BatchLens, INDEX, KEY, PROPERTY

@testset "PropertyBatchLens" begin
    l = PropertyBatchLens(:a, :b, :c)

    @test get((a=1, b=2, c=3), l) === (a=1, b=2, c=3)
    @test set((a=1, b=2, c=3), l, (a=10, b=20, c=30)) === (a=10, b=20, c=30)
end

@testset "KeyBatchLens" begin
    l = KeyBatchLens(:a, :b, :c)

    @test get((a=1, b=2, c=3), l) === (a=1, b=2, c=3)
    @test set((a=1, b=2, c=3), l, Dict(:a=>10, :b=>20, :c=>30)) === (a=10, b=20, c=30)
end

@testset "IndexBatchLens" begin
    l = IndexBatchLens(:a, :b, :c)

    @test get((a=1, b=2, c=3), l) === (1, 2, 3)
    @test get((b=2, a=1, c=3), l) === (1, 2, 3)
    @test get((b=2, c=3, a=1), l) === (1, 2, 3)
    @test set((a=1, b=2, c=3), l, (10, 20, 30)) === (a=10, b=20, c=30)
end

@testset "get" begin
    names = (:a, :b, :c)
    nt = (d=4, c=3, b=2, a=1)
    @testset "BatchLens{_, $objacc, _}" for (objacc, desired) in [
        (INDEX, (4, 3, 2))
        (KEY, (1, 2, 3))
        (PROPERTY, (1, 2, 3))
    ]
        lens = BatchLens{names, objacc, INDEX}()
        @test get(nt, lens) == desired
        if objacc === KEY
            @test get(Dict(pairs(nt)), lens) == desired
        end
    end
end

const Associative = Union{AbstractDict, NamedTuple}

eq(x, y) = x == y
eq(x::Associative, y::Associative) =
    sort(collect(keys(x))) == sort(collect(keys(y))) &&
    all(eq.(getindex.(Ref(x), keys(x)),
            getindex.(Ref(y), keys(x))))
eq(x::Tuple, y::NamedTuple) = x == Tuple(y)
eq(x::NamedTuple, y::Tuple) = Tuple(x) == y

@testset "laws" begin
    @testset "$obj" for (obj, objacc_list) in [
        ((a=1, b=2, c=3, d=4), [PROPERTY]),
        ((1, 2, 3, 4), [INDEX]),
        (Dict(:a=>1, :b=>2, :c=>3, :d=>4), [KEY]),
    ]
        @testset "BatchLens{_, $objacc, $valueacc}" for
                objacc in objacc_list,
                (valueacc, val_list) in [
                    (INDEX, [
                        (1, 2, 3),
                        (c=4, b=5, a=6),
                    ]),
                    (KEY, [
                        (a=1, b=2, c=3),
                        (c=4, b=5, a=6),
                        Dict(:a=>7, :b=>8, :c=>9),
                    ]),
                ]

            lens = BatchLens{(:a, :b, :c), objacc, valueacc}()

            @testset "You get what you set." begin
                @testset for val in val_list
                    @debug(
                        "eq(get(set(obj, lens, val), lens), val)",
                        obj,
                        lens,
                        lhs = get(set(obj, lens, val), lens),
                        rhs = val,
                    )
                    @test eq(get(set(obj, lens, val), lens), val)
                end
            end

            @testset "Setting what was already there changes nothing." begin
                @test set(obj, lens, get(obj, lens)) == obj
            end

            @testset "The last set wins." begin
                val1 = val_list[1]
                @testset for val2 in val_list
                    @test set(set(obj, lens, val1), lens, val2) ==
                        set(obj, lens, val2)
                end
            end
        end
    end
end

end  # module


text:
#--- Chain Rule

g(x)=(f3 *(f2(x)f1(x))
#=
chain rules give  #TODO: a loop on allJacobians

@nickhigham:
suppose you have a function g(of x) -
    a function of Composition
  (of other functions: f3, f2, f1)
  each f maps Vectors (of length nk -to-> nk+1)
1.so the
(1st) derivative of g = Jacobian Matrix
2.if Apply the chain Rule:
I'll get the expression (Jf):
the product of 3 Jacobians
Evaluated at the appropriate arguments

Question:
if I want to get the Jacobian (J)
  which order should I take
  J3 J2 J1 or J1, J2, J3 ?

  comes up in automatic differentiation:
  the difference between the first mode

Reverse mode (output) & forward (input )

so if i want to get the derivative of g

#-- testing Area:

for i in enumerate(n)

  if i <2
    return jf(i,)

=#
#Composition
Jg(x) =  Jf3(f2*(f1(x)))*Jf2(f1(x))*Jf1(x) ;

#---Todo implement Jacobian too

#=y Symbolics.substitute
https://discourse.julialang.org/t/how-to-evaluate-substitute-numerical-values-of-params-n-states-a-jacobian-in-modelingtoolkit-using-generate-jacobian-or-calculate-jacobian/71466

Solution
yewalenikhil65
Nov 2021
thanks @baggepinnen
i am doing currently following

states ,parameters
i == 1
x1map =

warning: ModelingToolkit uses ZygoteRules
plus next line it uses Symbolics too writer seems
can we do better ?
=#
#= Requires ModelingToolkit, Symbolics, Zygote (3 in 1)
pmap = nothing ;jac=nothing ;x1=0;ps =;
x₀map = states(odesys) .=> x1
pmap  = parameters(odesys) .=> ps
jac = substitute.( ModelingToolkit.calculate_jacobian(odesys), (Dict([x₀map;pmap]),) )
jac = Symbolics.value.(jac);
=#
#=the key
You only need function f(x)

Defining F as zeros(3) restricts its entries to be Float64s, whereas to use ForwardDiff
they need to be Duals.

Change it to F = similar(x) or F = zero.(x)
#https://discourse.julialang.org/t/jacobian-of-a-multivariate-function/21131/10
=#

#=@dpsanders

to get something of the correct type: =#

#--- chain rule gives
#=
function _f(x)
  F = similar(x) #zero.(x)#small perturbation

  F[1] = x[1]^2 + x[3] # function here
  F[2] = x[1] + x[2] #another function
  F[3] = x[2]^2 + x[3]^2 #third one too

  return F
end
=#
g(x) = f3(f2(f1(x)))   # fk R^n => R^nk+1;

Jg(x) = Jf3(f2(f1(x)))*Jf2(f1(x))*Jf1(x);

#=
Q. does Jg as (Jf3*Jf2)*Jf1 or Jf3*(Jf2*(Jf1))
me: i.e. is it forward or backward (check also nick-Higham juliacon 2018)

https://github.com/JuliaDiff/SparseDiffTools.jl
In addition, the following forms allow you to provide a gradient function g(dy,x) or dy=g(x) respectively:
=#
num_hesvecgrad!(dy,g,x,v,
                     cache2 = similar(v),
                     cache3 = similar(v))

num_hesvecgrad(g,x,v)

auto_hesvecgrad!(dy,g,x,v,
                     cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v),
                     cache3 = ForwardDiff.Dual{DeivVecTag}.(x, v))

auto_hesvecgrad(g,x,v);

#---
#The `numauto`and autonum methods both mix numerical and automatic differentiation, with the former almost always being more efficient and thus being recommended.

# Optionally, if you load Zygote.jl, the following numback and autoback methods are available and allow numerical/ForwardDiff over reverse mode automatic differentiation respectively, where the reverse-mode AD is provided by Zygote.jl. Currently these methods are not competitive against numauto, but as Zygote.jl gets optimized these will likely be the fastest.

#using Zygote # Required #i'm not sure # are there anytthing else?

numback_hesvec!(dy,f,x,v,
                     cache1 = similar(v),
                     cache2 = similar(v))

numback_hesvec(f,x,v)

# Currently errors! See https://github.com/FluxML/Zygote.jl/issues/241 #says about Zygote
autoback_hesvec!(dy,f,x,v,
                     cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v),
                     cache3 = ForwardDiff.Dual{DeivVecTag}.(x, v))

autoback_hesvec(f,x,v)
#=
Jv and Hv Operators
The following produce matrix-free operators which are used for calculating Jacobian-vector and Hessian-vector products where the differentiation takes place at the vector u:
=#
JacVec(f,x::AbstractArray;autodiff=true)
HesVec(f,x::AbstractArray;autodiff=true)
HesVecGrad(g,x::AbstractArray;autodiff=false)

#These all have the same interface, where J*v utilizes the out-of-place Jacobian-vector or Hessian-vector function, whereas mul!(res,J,v) utilizes the appropriate in-place versions. To update the location of differentiation in the operator, simply mutate the vector u: J.u .= ....

#=Chain Rule Modes of differentiation - Prof. Nick Higham

 1.Automatic Differentiation (AD)
 2.Forward (mode)
 3.Reverse mode

Provides a  references
G.Strang Linear Algebra & Learning from Data
Wellesley-Cambridge Press 2018
http://math.mit.edu/~gs/Learningfromdata

#Final inference:
no body has a clear idea about Automatic differentiation - and is still
a ripe  area of Research

"No body knows anything about anything"

=#

max_stars_repo_path: Tutorial/Nick Higham/composition.jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_path: Tutorial/Nick Higham/composition.jl

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_path: Tutorial/Nick Higham/composition.jl

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
__precompile__()
module ObjFileBase

export ObjectHandle, SectionRef, SymbolRef, debugsections

export printfield, printentry, printfield_with_color, deref,
    sectionaddress, sectionoffset, sectionsize, sectionname,
    load_strtab, readmeta, StrTab, symname, Sections, symbolvalue,
    isundef, symbolnum, Symbols, isBSS

import Base: read, seek, readbytes, position, show, showcompact, readuntil, skip

########## ObjFileBase.jl - Basic shared functionality for all object files ####
#
#   This package provides basic functionality that is shared among object files.
#   It was written with ELF/MachO/COFF in mind, but should be easily adaptable
#   to other object files as well.
#
#   The family of Object file implementations have a somewhat different focus than
#   a lot of other ones. While certainly usable to write a linker or compiler,
#   they are more intended for interactive exploration and verification and some
#   of the design choices reflect this. Nevertheless, this does not mean that
#   performance should be unnecessarily compromised.
#
# # # Basic Concepts
#
# # ObjectHandle
#
#   Provides an abstraction over a specific instance of the object file.
#   Currently many implementations put the IO object to be read from in this
#   handle, which has proven convenient for interactive usage.
#
# # SectionIterator, SymbolIterator, etc.
#
#   Provides support for iterating over Symbols, Sections, Relocations, etc.
#   These should have a reference to the original ObjectHandle and read the
#   desired information on demand. Rather than returning the read datastructure
#   directly, it should return in in an appropriate *Ref object.
#
# # SetionRef, SymbolRef
#
#   Convenience objects (mostly for interactive usage) that contain a reference
#   to the original object handle as well as any information that may be needed
#   to get further information about the Referenced object that is not
#   explicitly encoded in the data structure itself. Examples may include:
#
#       - Offset in the file (especially for variable length records
#            where only the header is part of the datastructure)
#       - Symbol/Section index
#       - Etc.
#
################################################################################

# # # Tools

macro mustimplement(sig)
    fname = sig.args[1]
    arg1 = sig.args[2]
    if isa(arg1,Expr)
        arg1 = arg1.args[1]
    end
    :($(esc(sig)) = error(typeof($(esc(arg1))),
                          " must implement ", $(Expr(:quote,fname))))
end

################################# Interfaces ###################################

# # # ObjectHandle

abstract ObjectHandle

# External facing interface
readmeta{T<:ObjectHandle}(io::IO, ::Type{T}) =
    error("$T must register a readmeta method")

@mustimplement readheader(oh::ObjectHandle)

# Querying general properties
@mustimplement endianness(oh::ObjectHandle)
@mustimplement intptr(oh::ObjectHandle)

# ObjectHandle IO interface
@mustimplement seek(oh::ObjectHandle, args...)
@mustimplement position(oh::ObjectHandle)

"""
Determine whether an object file is a relocatable (.o) object file.
"""
function isrelocatable
end

"""
Determine whether an object file is an executable object file.
"""
function isexecutable
end

"""
Turn a section name into the object-file specific naming convention.
"""
function mangle_sname
end

##
#  These are parameterized on the type of Object Handle. An imaginary Foo
#  fileformat might declare:
#
#   module Foo
#       immutable FooHandle; ... end
#       immutable Section <: ObjFileBase.Section{FooHandle}
#       ...; end
#   end
##

abstract Symbols{T<:ObjectHandle}
abstract Sections{T<:ObjectHandle}
abstract SectionRef{T<:ObjectHandle}
abstract Section{T<:ObjectHandle}
abstract Relocation{T<:ObjectHandle}
abstract RelocationRef{T<:ObjectHandle}

# The name of the section
@mustimplement sectionname(section::SectionRef)

# Retrieving the actual section datastructure
@mustimplement deref(section::SectionRef)
deref(section::Section) = section

# Retrieving the object handle
@mustimplement handle(section::SectionRef)

abstract SymbolRef{T<:ObjectHandle}
abstract SymtabEntry{T<:ObjectHandle}

"""
The size of the actual data contained in the section. This should exclude any
padding mandated by the file format e.g. due to alignment rules
"""
sectionsize(x::SectionRef) = sectionsize(deref(x))

"""
The address of the section in virtual memory.
"""
sectionaddress(x::SectionRef) = sectionaddress(deref(x))

"""
The offset of the section in the file.
"""
sectionoffset(x::SectionRef) = sectionoffset(deref(x))

handleT{T}(::Union{Type{SectionRef{T}}, Type{Section{T}}, Type{SymbolRef{T}},
    Type{SymtabEntry{T}}}) = T
Base.read(s::SectionRef, args...) = read(handle(s), args...)
Base.read{T}(s::SectionRef, ::Type{T}) = read(handle(s), T)::T
Base.position(s::SectionRef) = position(handle(s)) - sectionoffset(s)

abstract StrTab

function load_strtab
end
@mustimplement strtab_lookup(s::StrTab, offset)

# Symbol properties
function symname
end
function symbolvalue
end

symbolnum(x::SymbolRef) = symbolnum(deref(x))
isundef(x::SymbolRef) = isundef(deref(x))

# Section properties
isBSS(x::SectionRef) = isBSS(deref(x))

################################# Utilities ####################################

typealias SectionOrRef{T} Union{Section{T},SectionRef{T}}

seek{T<:ObjectHandle,S}(oh::T, section::Section{S}) =
    (@assert T <: S; seek(oh,sectionoffset(section)))

seek(section::SectionRef) = seek(handle(section), deref(section))
Base.eof(section::SectionRef) = eof(handle(section)) ||
    (position(handle(section)) >= sectionoffset(section) + sectionsize(section))

function readbytes{T<:ObjectHandle,S}(oh::T,sec::Section{S})
    @assert T <: S
    seek(oh, sec)
    readbytes(oh, sectionsize(sec))
end
readbytes(sec::SectionRef) = readbytes(handle(sec),deref(sec))
seek(x::SectionRef, off) = seek(handle(x), sectionoffset(x)+off)

typealias Maybe{T} Union{T,Void}

# # # Higher level debug info support
immutable DebugSections{T<:ObjectHandle, Sect}
    oh::T
    debug_abbrev::Maybe{Sect}
    debug_aranges::Maybe{Sect}
    debug_frame::Maybe{Sect}
    debug_info::Maybe{Sect}
    debug_line::Maybe{Sect}
    debug_loc::Maybe{Sect}
    debug_macinfo::Maybe{Sect}
    debug_pubnames::Maybe{Sect}
    debug_ranges::Maybe{Sect}
    debug_str::Maybe{Sect}
    debug_types::Maybe{Sect}
end
handle(dbgs::DebugSections) = dbgs.oh

function DebugSections{T}(oh::T; debug_abbrev = nothing, debug_aranges = nothing,
    debug_frame = nothing, debug_info = nothing, debug_line = nothing,
    debug_macinfo = nothing, debug_pubnames = nothing, debug_loc= nothing,
    debug_ranges = nothing, debug_str = nothing, debug_types = nothing)
    Sect = Union{map(typeof, [debug_abbrev, debug_aranges, debug_frame, debug_info,
        debug_line, debug_loc, debug_macinfo, debug_pubnames, debug_ranges,
        debug_str, debug_types])...}
    DebugSections{T,Sect}(oh, debug_abbrev, debug_aranges, debug_frame, debug_info,
        debug_line, debug_loc, debug_macinfo, debug_pubnames, debug_ranges,
        debug_str, debug_types)
end

function DebugSections{T<:ObjectHandle}(oh::T, sections::Dict)
    DebugSections(oh,
        debug_abbrev = get(sections, "debug_abbrev", nothing),
        debug_aranges = get(sections, "debug_aranges", nothing),
        debug_frame = get(sections, "debug_frame", nothing),
        debug_info = get(sections, "debug_info", nothing),
        debug_line = get(sections, "debug_line", nothing),
        debug_loc = get(sections, "debug_loc", nothing),
        debug_macinfo = get(sections, "debug_macinfo", nothing),
        debug_pubnames = get(sections, "debug_pubnames", nothing),
        debug_ranges = get(sections, "debug_ranges", nothing),
        debug_str = get(sections, "debug_str", nothing),
        debug_types = get(sections, "debug_types", nothing))
end

const DEBUG_SECTIONS = [
    "debug_abbrev",
    "debug_aranges",
    "debug_frame",
    "debug_info",
    "debug_line",
    "debug_loc",
    "debug_macinfo",
    "debug_pubnames",
    "debug_ranges",
    "debug_str",
    "debug_types"]

function show(io::IO, dsect::DebugSections)
    println(io, "Debug Sections for $(dsect.oh)")
    println(io,"========================= debug_abbrev =========================")
    println(io,dsect.debug_abbrev)
    println(io,"======================== debug_aranges =========================")
    println(io,dsect.debug_aranges)
    println(io,"========================= debug_frame ==========================")
    println(io,dsect.debug_frame)
    println(io,"========================= debug_info ===========================")
    println(io,dsect.debug_info)
    println(io,"========================= debug_line ===========================")
    println(io,dsect.debug_line)
    println(io,"========================== debug_loc ===========================")
    println(io,dsect.debug_loc)
    println(io,"======================== debug_macinfo =========================")
    println(io,dsect.debug_macinfo)
    println(io,"======================= debug_pubnames =========================")
    println(io,dsect.debug_pubnames)
    println(io,"======================== debug_ranges ==========================")
    println(io,dsect.debug_ranges)
    println(io,"=========================== debug_str ==========================")
    println(io,dsect.debug_str)
    println(io,"========================= debug_types ==========================")
    println(io,dsect.debug_types)
end

@mustimplement debugsections(oh::ObjectHandle)

export findindexbyname, findcubyname

function findindexbyname
end

function findcubyname
end

# Utils
# JIT Utils

function is_jit_section(s::Section)
    sectionaddress(s) > 0x100000
end
is_jit_section(s::SectionRef) = is_jit_section(deref(s))

@mustimplement replace_sections_from_memory(h::ObjectHandle, buffer)

# Printing utils
function printfield(io::IO,string,fieldlength; align = :right)
    (align == :right) && print(io," "^max(fieldlength-length(string),0))
    print(io,string)
    (align == :left) && print(io," "^max(fieldlength-length(string),0))
end
function printfield_with_color(color,io::IO,string,fieldlength; align = :right)
    (align == :right) && print(io," "^max(fieldlength-length(string),0))
    print_with_color(color,io,string)
    (align == :left) && print(io," "^max(fieldlength-length(string),0))
end
printentry(io::IO,header,values...) = (printfield(io,header,21);println(io," ",values...))

# User facing interfaces

immutable MagicMismatch;
    message
end

const ObjHandles = Vector{Type}()
function readmeta(io::IO)
    pos = position(io)
    for T in ObjHandles
        seek(io,pos)
        try
            return readmeta(io, T)
        catch e
            if !isa(e,MagicMismatch)
                rethrow(e)
            end
        end
    end
    error("""
        Object file is not any of $(join(ObjHandles, ", "))!
        To force one object file use readmeta(io,T).
        If the format you want is not listed, make sure
        the appropriate pacakge is loaded before calling
        this function.
        """)
end

readmeta(file::AbstractString) = readmeta(open(file,"r"))

read{T<:ObjectHandle}(oh::T, args...) = read(oh.io,args...)
read{H<:ObjectHandle,T}(oh::H, ::Type{T}) = read(oh.io, T)::T
readuntil{H<:ObjectHandle}(oh::H, x) = readuntil(oh.io, x)
skip{H<:ObjectHandle}(oh::H, x) = skip(oh.io, x)
# Others

function getSectionLoadAddress
end

immutable LOIByName
    addrs::Dict{Symbol, UInt64}
end
getSectionLoadAddress(LOI::LOIByName, x::Union{Symbol, AbstractString}) =
    LOI.addrs[symbol(x)]
getSectionLoadAddress(LOI::LOIByName, sec) =
    getSectionLoadAddress(LOI, bytestring(sectionname(sec)))

end # module


text:
module DataBags

export
    DataBag,
    contents,
    wrap

using Dates

"""

`AbstractDataBag{K,V,D}` is the super-type of data-bags.

"""
abstract type AbstractDataBag{K,V,D<:AbstractDict{K,V}} <: AbstractDict{K,V} end

struct DataBag{K,V,D<:AbstractDict{K,V}} <: AbstractDataBag{K,V,D}
    data::D
    DataBag{K,V,D}(data::D) where {K,V,D<:AbstractDict{K,V}} =
        new{K,V,D}(data)
end

# Generic outer constructors relying on `contents(Dict{K,V},...)` to build
# their contents.
DataBag(args...; kwds...) =
    wrap(DataBag, contents(Dict{Any,Any}, args...; kwds...))
DataBag{K}(args...; kwds...) where {K} =
    wrap(DataBag, contents(Dict{K,Any}, args...; kwds...))
DataBag{K,V}(args...; kwds...) where {K,V} =
    wrap(DataBag, contents(Dict{K,V}, args...; kwds...))

# Extends the `contents` methods to benefit from the API of `AbstractDataBag`.
@inline contents(A::DataBag) = Base.getfield(A, :data)

"""
    DataBags.contents(A)

yields the contents associated with an instance `A` of a sub-type of
`DataBags.AbstractDataBag` or of `AbstractDict`.

This method should be specialized by types derived from
`DataBags.AbstractDataBag`, this is the most simple way to inherit of the
common behavior implemented by this abstract type.

The `contents` method is also meant to be called by data-bag constructors to
create a new dictionary out of their arguments.  For that usage, the syntax is:

    content(Dict{K,V}, args...; kwds...) -> Dict

which yields a new dictionary built out of arguments `args` and keywords `kwds`
and accounting for type constraints set by `K` and `V` for respectively the
keys and values of the returned dictionary.  Types `K` and/or `V` can be ``Any`
if no type constraints are imposed.

""" contents

@noinline contents(::T) where {T<:AbstractDataBag} =
    error(string("method DataBags.contents has not been specialized for ", T))
@inline contents(A::AbstractDict) = A

# Methods for initial contents set from an existing dictionary data
# (e.g. another data-bag).  We create a new empty dictionary with proper
# key and value types and fill it with the contents of the argument.
contents(::Type{Dict{Any,Any}}, A::AbstractDict{K}) where {K} =
    merge!(Dict{K,Any}(), contents(A))
contents(::Type{Dict{K,Any}}, A::AbstractDict) where {K} =
    merge!(Dict{K,Any}(), contents(A))
contents(::Type{Dict{Any,V}}, A::AbstractDict{K}) where {K,V} =
    merge!(Dict{K,V}(), contents(A))
contents(::Type{Dict{K,V}}, A::AbstractDict) where {K,V} =
    merge!(Dict{K,V}(), contents(A))

# Methods for initial contents specified by keywords.
contents(::Type{Dict{Any,Any}}; kwds...) = Dict{Symbol,Any}(kwds...)
contents(::Type{Dict{Symbol,Any}}; kwds...) = Dict{Symbol,Any}(kwds...)
contents(::Type{Dict{Any,V}}; kwds...) where {V} = Dict{Symbol,V}(kwds...)
contents(::Type{Dict{Symbol,V}}; kwds...) where {V} = Dict{Symbol,V}(kwds...)

# Methods for initial contents specified as key-value pairs.  Unless key or
# value types are explictely specified, we want to have a dictionary whose
# key type is the most specific (for efficiency) and having the least
# specific value type (for flexibility).  If type constraint or type of
# pairs are insufficient, we eventually first make a dictionary and
# convert it so as to use the least specific value type.  Conversion is
# no-op if the dictionary is already of the correct type.  Since the key
# type must be known, we need a helper function of the initial dictionary
# must be known
contents(::Type{Dict{K,V}}, args::Pair...) where {K,V} = Dict{K,V}(args...)
contents(::Type{Dict{K,Any}}, args::Pair...) where {K} = Dict{K,Any}(args...)
contents(::Type{Dict{Any,V}}, args::Pair{K}...) where {K,V} = Dict{K,V}(args...)
contents(::Type{Dict{Any,V}}, args::Pair{<:AbstractString}...) where {V} =
    Dict{String,V}(args...)
contents(::Type{Dict{Any,V}}, args::Pair...) where {V} =
    withspecificvaluetype(V, Dict(args...))
contents(::Type{Dict{Any,Any}}, args::Pair{K}...) where {K} =
    Dict{K,Any}(args...)
contents(::Type{Dict{Any,Any}}, args::Pair{<:AbstractString}...) =
    Dict{String,Any}(args...)
contents(::Type{Dict{Any,Any}}, args::Pair...) =
    withspecificvaluetype(Any, Dict(args...))

withspecificvaluetype(::Type{V}, data::Dict{<:Any,V}) where {V} = data
withspecificvaluetype(::Type{V}, data::Dict{K}) where {K,V} = Dict{K,V}(data)

"""
    wrap(T::Type{<:AbstractDataBag}, arg) -> obj::T

wraps argument `arg` in a data-bag of type `T`.  The returned object shares
its contents with `arg`.

This method can be used to build a data-bag from an existing dictionary
without duplicating the dictionary.  As an example:

    dict = Dict{Symbol,Any}(:a => 1, :foo => "bar")
    cont = wrap(DataBag, dict)
    cont.b = 33
    dict["b"] == 33 # yields `true`

If the statement `wrap(DataBag, dict)` is replaced by `DataBag(dict)`, the
result of the last statement is `false`.

"""
wrap(::Type{DataBag}, data::D) where {K,V,D<:AbstractDict{K,V}} =
    DataBag{K,V,D}(data)
wrap(::Type{DataBag{K}}, data::D) where {K,V,D<:AbstractDict{K,V}} =
    wrap(DataBag, data)
wrap(::Type{DataBag{K,V}}, data::D) where {K,V,D<:AbstractDict{K,V}} =
    wrap(DataBag, data)
wrap(::Type{DataBag{K,V,D}}, data::D) where {K,V,D<:AbstractDict{K,V}} =
    wrap(DataBag, data)

"""
    DataBags.@newtype T

creates a new data-bag type `T` which is a sub-type of
`DataBags.AbstractDataBag` with symbolic keys and values of `Any` type.

"""
macro newtype(sym)
    isa(sym, Symbol) || throw(ArgumentError("argument must be a symbol"))
    T = esc(sym)
    quote
        struct $T <: DataBags.AbstractDataBag{Symbol,Any,Dict{Symbol,Any}}
            data::Dict{Symbol,Any}
            $T(args...; kwds...) =
                new(DataBags.contents(Dict{Symbol,Any}, args...; kwds...))
        end
        @inline DataBags.contents(A::$T) = Base.getfield(A, :data)
    end
end

"""
    DataBags.propertyname(T, sym)

converts symbol `sym` to a suitable key for data-bag of type `T` (a sub-type of
`DataBags.AbstractDataBag`), throwing an error if this conversion is not
supported.

"""
propertyname(::Type{<:AbstractDataBag{Symbol}}, sym::Symbol) = sym
propertyname(::Type{<:AbstractDataBag{String}}, sym::Symbol) = String(sym)
@noinline propertyname(::Type{T}, sym::Symbol) where {K,T<:AbstractDataBag{K}} =
    error(string("Converting symbolic key to type `", K, "` is not ",
                 "implemented. As a result, syntax `obj.key` is not ",
                 "supported for objects of type `", T, "`."))

# Extend methods so that syntax `obj.field` can be used.
@inline Base.getproperty(A::T, sym::Symbol) where {T<:AbstractDataBag} =
    getindex(contents(A), propertyname(T, sym))
@inline Base.setproperty!(A::T, sym::Symbol, val) where {T<:AbstractDataBag} =
    setindex!(contents(A), val, propertyname(T, sym))

# FIXME: Should return `Tuple(keys(contents(A)))` to conform to the doc. of
#        `propertynames` but this is slower and for most purposes, an iterable
#        is usually needed.
Base.propertynames(A::AbstractDataBag, private::Bool=false) = keys(A)

Base.convert(::Type{T}, A::T) where {T<:AbstractDataBag} = A
Base.convert(::Type{T}, A::AbstractDataBag) where {T<:AbstractDataBag} =
    T(A)

Base.keytype(::T) where {T<:AbstractDataBag} = keytype(T)
Base.keytype(::Type{<:AbstractDataBag{K,V}}) where {K,V} = K

Base.valtype(::T) where {T<:AbstractDataBag} = valtype(T)
Base.valtype(::Type{<:AbstractDataBag{K,V}}) where {K,V} = V

Base.length(A::AbstractDataBag) = length(contents(A))

Base.iterate(A::AbstractDataBag) = iterate(contents(A))
Base.iterate(A::AbstractDataBag, state) = iterate(contents(A), state)

Base.haskey(A::AbstractDataBag, key) = haskey(contents(A), key)
Base.keys(A::AbstractDataBag) = keys(contents(A))
Base.values(A::AbstractDataBag) = values(contents(A))
Base.getkey(A::AbstractDataBag, key, def) = getkey(contents(A), key, def)
Base.get(A::AbstractDataBag, key, def) = get(contents(A), key, def)
Base.get!(A::AbstractDataBag, key, def) = get!(contents(A), key, def)
Base.getindex(A::AbstractDataBag, key) = getindex(contents(A), key)
Base.setindex!(A::AbstractDataBag, val, key) =
    (setindex!(contents(A), val, key); return A)
Base.delete!(A::AbstractDataBag, key) = begin
    delete!(contents(A), key)
    return A
end
Base.pop!(A::AbstractDataBag, key) = pop!(contents(A), key)
Base.pop!(A::AbstractDataBag, key, def) = pop!(contents(A), key, def)
Base.pairs(A::AbstractDataBag) = pairs(contents(A))

# Standard methods:
#
# - merge() and merge!() do not need to be specialized.
#
# - Specializing empty() and empty!() is sufficient for copy() and copy!()
#   to work.

Base.empty(A::DataBag{K,V}) where {K,V} =
    wrap(DataBag, empty(contents(A), K, V))
Base.empty(A::DataBag{<:Any,V}, ::Type{K}) where {K,V} =
    wrap(DataBag, empty(contents(A), K, V))
Base.empty(A::DataBag, ::Type{K}, ::Type{V}) where {K,V} =
    wrap(DataBag, empty(contents(A), K, V))

Base.empty!(A::AbstractDataBag) = begin
    empty!(contents(A))
    return A
end

function Base.show(io::IO, ::MIME"text/plain",
                   A::AbstractDataBag{K,V}) where {K,V}
    n = length(A)
    summary(io, A)
    println(io, ":")
    lstkeys = sort(collect(keys(A)))
    strkeys = repr.(lstkeys)
    maxlen = maximum(length, strkeys)
    spaces = [repeat(" ", r) for r in 0:maxlen]
    for i in 1:n
        key = lstkeys[i]
        str = strkeys[i]
        print(io, "  ", str, spaces[maxlen + 1 - length(str)], " => ")
        _showval(io, A[key])
        if (n -= 1) > 0
            println(io)
        end
    end
end

_showval(io::IO, val::Number) = print(io, val)
_showval(io::IO, val::AbstractString) = print(io, repr(val))
_showval(io::IO, val::AbstractRange) = print(io, repr(val))
_showval(io::IO, val::T) where {T} = summary(io, val)
_showval(io::IO, val::DateTime) = print(io, "DateTime(\"", repr(val), "\")")

end # module