text:
function volume_fractions!(fracs, topo, mesh, C)
    map!(fracs, 1:nelements(mesh)) do i
        el = element(mesh, i)
        T = typeof(el)
        inds = Tuple(indices(element(topo, i)))
        vs = Tuple(vertices(el))
        Cs = C[SVector(inds)]
        return positive_measure_ratio(T, inds, vs, Cs)
    end
end

function positive_measure_ratio(::Type{<:Triangle}, inds, vs, C)
    a1, a2 = _positive_measures(Triangle, inds, vs, C)
    return a1 / a2
end
function _positive_measures(::Type{<:Triangle}, inds, vs, C)
    nneg = 0
    npos = 0
    all_pos = true
    all_neg = true
    for ind in (1, 2, 3)
        if C[ind] < 0
            all_pos = false
            nneg += 1
        end
        if C[ind] > 0
            all_neg = false
            npos += 1
        end
    end
    a1 = Meshes.measure(Triangle(vs...)) |> abs
    if all_neg
        return zero(eltype(C)), a1
    elseif all_pos
        return convert(eltype(C), a1), a1
    end

    ind0 = -1
    for j in (1, 2, 3)
        if nneg == 1 && C[j] < 0
            ind0 = j
            break
        elseif nneg == 2 && C[j] > 0
            ind0 = j
            break
        end
    end
    @assert ind0 != -1
    others = setdiff((1, 2, 3), ind0)

    δC1 = C[ind0] - C[others[1]]
    δx1 = vs[ind0] - vs[others[1]]
    ratio1 = C[ind0] / δC1 * norm(δx1)
    x1 = vs[ind0] - ratio1 * δx1

    δC2 = C[ind0] - C[others[2]]
    δx2 = vs[ind0] - vs[others[2]]
    ratio2 = C[ind0] / δC2 * norm(δx2)
    x2 = vs[ind0] - ratio2 * δx2

    TP = typeof(x1)
    a2 = Meshes.measure(Triangle(x1, x2, convert(TP, vs[ind0]))) |> abs

    if nneg == 1
        return (a1 - a2), a1
    else
        return a2, a1
    end
end

function positive_measure_ratio(::Type{<:Quadrangle}, inds, vs, C)
    a1, a2 = _positive_measures(Quadrangle, inds, vs, C)
    return a1 / a2
end
function _positive_measures(::Type{<:Quadrangle}, inds, vs, C)
    inds1 = @SVector [1, 2, 3]
    a1, b1 = _positive_measures(Triangle, inds[inds1], vs[inds1], C[inds1])
    inds2 = @SVector [1, 3, 4]
    a2, b2 = _positive_measures(Triangle, inds[inds2], vs[inds2], C[inds2])
    return (a1 + a2), (b1 + b2)
end

function positive_measure_ratio(::Type{<:Tetrahedron}, inds, vs, C)
    a1, a2 = _positive_measures(Tetrahedron, inds, vs, C)
    return a1 / a2
end
function _positive_measures(::Type{<:Tetrahedron}, inds, vs, C)
    nneg = 0
    npos = 0
    all_pos = true
    all_neg = true
    for ind in (1, 2, 3, 4)
        if C[ind] < 0
            all_pos = false
            nneg += 1
        end
        if C[ind] > 0
            all_neg = false
            npos += 1
        end
    end
    a1 = Meshes.measure(Tetrahedron(vs...)) |> abs
    if all_neg
        return zero(eltype(C)), a1
    elseif all_pos
        return convert(eltype(C), a1), a1
    end

    neg_inds = Int[]
    for j in (1, 2, 3, 4)
        if C[j] < 0
            push!(neg_inds, j)
        end
    end
    if nneg == 1
        others = setdiff((1, 2, 3, 4), neg_inds)
        ind0 = neg_inds[1]

        δC1 = C[ind0] - C[others[1]]
        δx1 = vs[ind0] - vs[others[1]]
        ratio1 = C[ind0] / δC1 * norm(δx1)
        x1 = vs[ind0] - ratio1 * δx1

        δC2 = C[ind0] - C[others[2]]
        δx2 = vs[ind0] - vs[others[2]]
        ratio2 = C[ind0] / δC2 * norm(δx2)
        x2 = vs[ind0] - ratio2 * δx2

        δC3 = C[ind0] - C[others[3]]
        δx3 = vs[ind0] - vs[others[3]]
        ratio3 = C[ind0] / δC3 * norm(δx3)
        x3 = vs[ind0] - ratio3 * δx3

        TP = typeof(x1)
        a2 = Meshes.measure(Tetrahedron(x1, x2, x3, convert(TP, vs[ind0]))) |> abs

        return (a1 - a2), a1
    elseif length(neg_inds) == 2
        others = setdiff((1, 2, 3, 4), neg_inds)
        ind0, ind1 = neg_inds

        δC1 = C[ind0] - C[others[1]]
        δx1 = vs[ind0] - vs[others[1]]
        ratio1 = C[ind0] / δC1 * norm(δx1)
        x1 = vs[ind0] - ratio1 * δx1

        δC2 = C[ind0] - C[others[2]]
        δx2 = vs[ind0] - vs[others[2]]
        ratio2 = C[ind0] / δC2 * norm(δx2)
        x2 = vs[ind0] - ratio2 * δx2

        δC3 = C[ind1] - C[others[1]]
        δx3 = vs[ind1] - vs[others[1]]
        ratio3 = C[ind1] / δC3 * norm(δx3)
        x3 = vs[ind1] - ratio3 * δx3

        δC4 = C[ind1] - C[others[2]]
        δx4 = vs[ind1] - vs[others[2]]
        ratio4 = C[ind1] / δC4 * norm(δx4)
        x4 = vs[ind1] - ratio4 * δx4

        TP = typeof(x1)
        a2 = Meshes.measure(Wedge(x1, x2, convert(TP, vs[ind0]), x3, x4, convert(TP, vs[ind1]))) |> abs

        return (a1 - a2), a1
    elseif length(neg_inds) == 3
        pos_inds = setdiff((1, 2, 3, 4), neg_inds)
        ind0 = pos_inds[1]
        others = neg_inds

        δC1 = C[ind0] - C[others[1]]
        δx1 = vs[ind0] - vs[others[1]]
        ratio1 = C[ind0] / δC1 * norm(δx1)
        x1 = vs[ind0] - ratio1 * δx1

        δC2 = C[ind0] - C[others[2]]
        δx2 = vs[ind0] - vs[others[2]]
        ratio2 = C[ind0] / δC2 * norm(δx2)
        x2 = vs[ind0] - ratio2 * δx2

        δC3 = C[ind0] - C[others[3]]
        δx3 = vs[ind0] - vs[others[3]]
        ratio3 = C[ind0] / δC3 * norm(δx3)
        x3 = vs[ind0] - ratio3 * δx3

        TP = typeof(x1)
        a2 = Meshes.measure(Tetrahedron(x1, x2, x3, convert(TP, vs[ind0]))) |> abs

        return a2, a1
    else
        throw("Unreachable.")
    end
end

function positive_measure_ratio(::Type{<:Hexahedron}, inds, vs, C)
    a1, a2 = _positive_measures(Hexahedron, inds, vs, C)
    return a1 / a2
end
function _positive_measures(::Type{<:Hexahedron}, inds, vs, C)
    _inds = SVector(1, 2, 4, 5)
    a1, b1 = _positive_measures(Tetrahedron, inds[_inds], vs[_inds], C[_inds])
    _inds = SVector(2, 7, 4, 5)
    a2, b2 = _positive_measures(Tetrahedron, inds[_inds], vs[_inds], C[_inds])
    _inds = SVector(4, 5, 7, 8)
    a3, b3 = _positive_measures(Tetrahedron, inds[_inds], vs[_inds], C[_inds])
    _inds = SVector(2, 6, 7, 5)
    a4, b4 = _positive_measures(Tetrahedron, inds[_inds], vs[_inds], C[_inds])
    _inds = SVector(4, 7, 2, 3)
    a5, b5 = _positive_measures(Tetrahedron, inds[_inds], vs[_inds], C[_inds])
    return (a1 + a2 + a3 + a4 + a5), (b1 + b2 + b3 + b4 + b5)
end


text:
module ArrayInterfaceOffsetArrays

using ArrayInterface
using OffsetArrays
using Static

relative_offsets(r::OffsetArrays.IdOffsetRange) = (getfield(r, :offset),)
relative_offsets(A::OffsetArrays.OffsetArray) = getfield(A, :offsets)
function relative_offsets(A::OffsetArrays.OffsetArray, ::StaticInt{dim}) where {dim}
    if dim > ndims(A)
        return static(0)
    else
        return getfield(relative_offsets(A), dim)
    end
end
function relative_offsets(A::OffsetArrays.OffsetArray, dim::Int)
    if dim > ndims(A)
        return 0
    else
        return getfield(relative_offsets(A), dim)
    end
end
ArrayInterface.parent_type(::Type{<:OffsetArrays.OffsetArray{T,N,A}}) where {T,N,A} = A
function _offset_axis_type(::Type{T}, dim::StaticInt{D}) where {T,D}
    OffsetArrays.IdOffsetRange{Int,ArrayInterface.axes_types(T, dim)}
end
function ArrayInterface.axes_types(::Type{T}) where {T<:OffsetArrays.OffsetArray}
    Static.eachop_tuple(
        _offset_axis_type,
        ntuple(static, StaticInt(ndims(T))),
        ArrayInterface.parent_type(T)
    )
end
ArrayInterface.strides(A::OffsetArray) = ArrayInterface.strides(parent(A))
function ArrayInterface.known_offsets(::Type{A}) where {A<:OffsetArrays.OffsetArray}
    ntuple(identity -> nothing, Val(ndims(A)))
end
function ArrayInterface.offsets(A::OffsetArrays.OffsetArray)
    map(+, ArrayInterface.offsets(parent(A)), relative_offsets(A))
end
@inline function ArrayInterface.offsets(A::OffsetArrays.OffsetArray, dim)
    d = ArrayInterface.to_dims(A, dim)
    ArrayInterface.offsets(parent(A), d) + relative_offsets(A, d)
end
@inline function ArrayInterface.axes(A::OffsetArrays.OffsetArray)
    map(OffsetArrays.IdOffsetRange, ArrayInterface.axes(parent(A)), relative_offsets(A))
end
@inline function ArrayInterface.axes(A::OffsetArrays.OffsetArray, dim)
    d = ArrayInterface.to_dims(A, dim)
    OffsetArrays.IdOffsetRange(ArrayInterface.axes(parent(A), d), relative_offsets(A, d))
end
function ArrayInterface.stride_rank(T::Type{<:OffsetArray})
  ArrayInterface.stride_rank(ArrayInterface.parent_type(T))
end
function ArrayInterface.dense_dims(T::Type{<:OffsetArray})
    ArrayInterface.dense_dims(ArrayInterface.parent_type(T))
end
function ArrayInterface.contiguous_axis(T::Type{<:OffsetArray})
  ArrayInterface.contiguous_axis(ArrayInterface.parent_type(T))
end

end # module

ext: jl


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# CSG

"""
    valid_intervals(ts::Vector)

Return a `Vector` of `Interval`s describing what ranges of the `t` parameter are internal to the shape, given its `ts` hit parameters.

For the function to work correctly `ts` must be of even length.
"""
function valid_intervals(ts::Vector)
    iseven(length(ts)) || throw(ArgumentError("given vector must have even length: got $ts"))
    sorted_ts = sort(ts)
    [Interval{Open, Open}(t1, t2) for (t1, t2) ∈ zip(sorted_ts[begin:2:end], sorted_ts[begin+1:2:end])]
end

#######
# CSG

"""
    Rule

Enum type representing the hit point selection of a [`CSG`](@ref).

# Instances

- `UniteRule`: indicates that every hit point is valid
- `IntersectRule`: indicates that only hit points located inside of other shapes are valid
- `DiffRule`: indicates that only hit points outside of the `lbranch` and inside the `rbranch` are valid
- `FuseRule`: indicates that every hit point outside of other shapes is valid

See also: [`CSG`](@ref).
"""
@enum Rule begin
    UniteRule
    IntersectRule
    DiffRule
    FuseRule
end

"""
    CSG{R} <: CompositeShape

A [`Shape`](@ref) representing a Constructive Solid Geometry tree.

The behavior of the CSG tree is determined by the [`Rule`](@ref) `R`.

# Members

- `rbranch::Shape`: represents the right branch of the tree
- `lbranch::Shape`: represents the left branch of the tree
- `transformation::Transformation`: represents the [`Transformation`](@ref) of the whole composite shape

# External references

- Constructive Solid Geometry: https://en.wikipedia.org/wiki/Constructive_solid_geometry

"""
struct CSG{R} <: CompositeShape
    rbranch::Shape
    lbranch::Shape
    transformation::Transformation
    """
        CSG{R::Rule}(rbranch::Shape, lbranch::Shape)

    Constrcts an instance of a CSG tree with the given [`Rule`](@ref) and branches.
    """
    function CSG{R}(rbranch::Shape, lbranch::Shape; transformation::Transformation = Transformation()) where {R}
        R::Rule
        new{R}(rbranch, lbranch, transformation)
    end
end

"""
    UnionCSG

Alias for `CSG{UniteRule}`.

See also: [`CSG`](@ref), [`Rule`](@ref).
"""
const UnionCSG = CSG{UniteRule}

"""
    IntersectionCSG

Alias for `CSG{IntersectRule}`.

See also: [`CSG`](@ref), [`Rule`](@ref).
"""
const IntersectionCSG = CSG{IntersectRule}

"""
    DiffCSG

Alias for `CSG{DiffRule}`.

See also: [`CSG`](@ref), [`Rule`](@ref).
"""
const DiffCSG = CSG{DiffRule}

"""
    FusionCSG

Alias for `CSG{FuseRule}`.

See also: [`CSG`](@ref), [`Rule`](@ref).
"""
const FusionCSG = CSG{FuseRule}

"""
    union(s1::Shape, s2::Shape; transformation::Transformation = Transformation())

Construct a [`UnionCSG`](@ref) with the given shapes as `rbranch` and `lbranch` repectively.
"""
function Base.union(s1::Shape, s2::Shape; transformation::Transformation = Transformation())
    UnionCSG(s1, s2, transformation = transformation)
end

"""
    union(s::Shape, ss::Shape...; transformation::Transformation = Transformation())

Construct a [`UnionCSG`](@ref) binary tree, by recursively calling [`union`](@ref)`(::Shape, ::Shape)`.
"""
function Base.union(s::Shape, ss::Shape...; transformation::Transformation = Transformation())
    union(union(s, ss[begin:end ÷ 2]...), union(ss[end ÷ 2 + 1:end]...); transformation = transformation)
end

function Base.union(s1::Shape, s2::Shape, s3::Shape; transformation::Transformation = Transformation())
    union(union(s1, s2), s3; transformation = transformation)
end

"""
    intersect(s1::Shape, s2::Shape; transformation::Transformation = Transformation())

Construct a [`IntersectionCSG`](@ref) with the given shapes as `rbranch` and `lbranch` repectively.
"""
function Base.intersect(s1::Shape, s2::Shape; transformation::Transformation = Transformation())
    IntersectionCSG(s1, s2; transformation = transformation)
end

"""
    intersect(s::Shape, ss::Shape...; transformation::Transformation = Transformation())

Construct a [`IntersectionCSG`](@ref) binary tree, by recursively calling [`intersect`](@ref)`(::Shape, ::Shape)`.
"""
function Base.intersect(s::Shape, ss::Shape...; transformation::Transformation = Transformation())
    intersect(intersect(s, ss[begin:end ÷ 2]...), intersect(ss[(end ÷ 2 + 1):end]...); transformation = transformation)
end

function Base.intersect(s1::Shape, s2::Shape, s3::Shape; transformation::Transformation = Transformation())
    intersect(intersect(s1, s2), s3; transformation = transformation)
end

"""
    setdiff(s1::Shape, s2::Shape); transformation::Transformation = Transformation())

Construct a [`DiffCSG`](@ref) with the given shapes as `rbranch` and `lbranch` repectively.
"""
function Base.setdiff(s1::Shape, s2::Shape; transformation::Transformation = Transformation())
    DiffCSG(s1, s2; transformation = transformation)
end

"""
    setdiff(s::Shape, ss::Shape...); transformation::Transformation = Transformation())

Construct a [`DiffCSG`](@ref) between `s` and [`fuse`](@ref)`(ss...)`.
"""
function Base.setdiff(s::Shape, ss::Shape...; transformation::Transformation = Transformation())
    setdiff(s, fuse(ss...), transformation = transformation)
end

"""
    fuse(s1::Shape, s2::Shape); transformation::Transformation = Transformation())

Construct a [`FusionCSG`](@ref) with the given shapes as `rbranch` and `lbranch` repectively.
"""
function fuse(s1::Shape, s2::Shape; transformation::Transformation = Transformation())
    FusionCSG(s1, s2; transformation = transformation)
end

"""
    fuse(s::Shape, ss::Shape...); transformation::Transformation = Transformation())

Construct a [`FusionCSG`](@ref) binary tree, by recursively calling [`intersect`](@ref)`(::Shape, ::Shape)`.
"""
function fuse(s::Shape, ss::Shape...; transformation::Transformation = Transformation())
    fuse(fuse(s, ss[begin:end ÷ 2]...), fuse(ss[end ÷ 2 + 1:end]...); transformation = transformation)
end

function fuse(s1::Shape, s2::Shape, s3::Shape; transformation::Transformation = Transformation())
    fuse(fuse(s1, s2), s3; transformation = transformation)
end

###########
# UnionCSG

function ray_intersection(ray::Ray, csg::UnionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_hit = ray_intersection(inv_ray, csg.rbranch)
    l_hit = ray_intersection(inv_ray, csg.lbranch)
    isnothing(r_hit) && return l_hit
    isnothing(l_hit) && return r_hit
    min(r_hit, l_hit)
end

function all_ray_intersections(ray::Ray, csg::UnionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_hits = all_ray_intersections(inv_ray, csg.rbranch)
    append!(r_hits , filter(hit -> hit ∉ r_hits, all_ray_intersections(inv_ray, csg.lbranch)))
end

function quick_ray_intersection(ray::Ray, csg::UnionCSG)
    inv_ray = inv(csg.transformation) * ray
    quick_ray_intersection(inv_ray, csg.rbranch) || quick_ray_intersection(inv_ray, csg.lbranch)
end

function get_all_ts(csg::UnionCSG, ray::Ray)
    inv_ray = inv(csg.transformation) * ray
    append!(get_all_ts(csg.rbranch, inv_ray), get_all_ts(csg.lbranch, inv_ray))
end

##################
# IntersectionCSG

function ray_intersection(ray::Ray, csg::IntersectionCSG)
    hits = filter(hit -> ray.tmin < hit.t < ray.tmax, all_ray_intersections(ray, csg))
    isempty(hits) && return nothing
    minimum(hits)
end

function all_ray_intersections(ray::Ray, csg::IntersectionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_hits = all_ray_intersections(inv_ray, csg.rbranch)
    isempty(r_hits) && return Vector{HitRecord}()
    l_hits = all_ray_intersections(inv_ray, csg.lbranch)
    isempty(l_hits) && return Vector{HitRecord}()
    r_intervals = valid_intervals(r_hits)
    l_intervals = valid_intervals(l_hits)
    r_filter = filter(hit -> any(Ref(hit) .∈ l_intervals), r_hits)
    l_filter = filter(hit -> any(Ref(hit) .∈ r_intervals) && hit ∉ r_filter, l_hits)
    append!(r_filter, l_filter)
end

function quick_ray_intersection(ray::Ray, csg::IntersectionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    isempty(r_ts) && return false
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(l_ts) && return false
    r_intervals = valid_intervals(r_hits)
    l_intervals = valid_intervals(l_hits)
    any(t -> any(t .∈ l_intervals) && ray.tmin < t < ray.tmax, r_ts) ||
    any(t -> any(t .∈ r_intervals) && ray.tmin < t < ray.tmax, l_ts)
end

function get_all_ts(csg::IntersectionCSG, ray::Ray)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    isempty(r_ts) && return Vector{Float32}()
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(l_ts) && return Vector{Float32}()
    r_intervals = valid_intervals(r_ts)
    l_intervals = valid_intervals(l_ts)
    r_filter = filter(t -> any(t .∈ l_intervals), r_ts)
    l_filter = filter(t -> any(t .∈ r_intervals) && t ∉ r_filter, l_ts)
    # @assert (length(r_filter) + length(l_filter) != 1) "Only one intersection for $(typeof(csg.rbranch)): $r_ts + $l_ts between $l_min and $l_max"
    append!(r_filter, l_filter)
end

##########
# DiffCSG

function ray_intersection(ray::Ray, csg::DiffCSG)
    hits = filter(hit -> ray.tmin < hit.t < ray.tmax, all_ray_intersections(ray, csg))
    isempty(hits) && return nothing
    minimum(hits)
end

function all_ray_intersections(ray::Ray, csg::DiffCSG)
    inv_ray = inv(csg.transformation) * ray
    r_hits = all_ray_intersections(inv_ray, csg.rbranch)
    isempty(r_hits) && return Vector{HitRecord}()
    l_hits = all_ray_intersections(inv_ray, csg.lbranch)
    isempty(l_hits) && return r_hits
    r_intervals = valid_intervals(r_hits)
    l_intervals = valid_intervals(l_hits)
    r_filter = filter(hit -> all(Ref(hit) .∉ l_intervals), r_hits)
    l_filter = filter(hit -> any(Ref(hit) .∈ r_intervals) && hit ∉ r_filter, l_hits)
    append!(r_filter, l_filter)
end

function quick_ray_intersection(ray::Ray, csg::DiffCSG)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    isempty(r_ts) && return false
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(l_ts) && return any(t -> ray.tmin < t < ray.tmax, r_ts)

    r_intervals = valid_intervals(r_ts)
    l_intervals = valid_intervals(l_ts)
    any(t -> all(t .∉ l_intervals) && ray.tmin < t < ray.tmax, r_ts) ||
    any(t -> any(t .∈ r_intervals) && ray.tmin < t < ray.tmax, l_ts)
end

function get_all_ts(csg::DiffCSG, ray::Ray)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    isempty(r_ts) && return Vector{Float32}()
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(l_ts) && return r_ts
    # r_min, r_max = extrema(r_hits)
    # l_min, l_max = extrema(l_hits)
    r_intervals = valid_intervals(r_ts)
    l_intervals = valid_intervals(l_ts)
    r_filter = filter(t -> any(t .∉ l_intervals), r_ts)
    l_filter = filter(t -> any(t .∈ r_intervals) && t ∉ r_filter, l_ts)
    append!(r_filter, l_filter)
end

############
# FusionCSG

function ray_intersection(ray::Ray, csg::FusionCSG)
    hits = filter(hit -> ray.tmin < hit.t < ray.tmax, all_ray_intersections(ray, csg))
    isempty(hits) && return nothing
    minimum(hits)
end

function all_ray_intersections(ray::Ray, csg::FusionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_hits = all_ray_intersections(inv_ray, csg.rbranch)
    l_hits = all_ray_intersections(inv_ray, csg.lbranch)
    isempty(r_hits) && return l_hits
    isempty(l_hits) && return r_hits
    r_intervals = valid_intervals(r_hits)
    l_intervals = valid_intervals(l_hits)
    r_filter = filter(hit -> all(Ref(hit) .∉ l_intervals), r_hits)
    l_filter = filter(hit -> all(Ref(hit) .∉ r_intervals) && hit ∉ r_filter, l_hits)
    isempty(r_filter) && return l_filter
    isempty(l_filter) && return r_filter
    append!(r_filter, l_filter)
end

function quick_ray_intersection(ray::Ray, csg::FusionCSG)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(r_ts) && return !isempty(l_ts)
    isempty(l_ts) && return true
    r_intervals = valid_intervals(r_ts)
    l_intervals = valid_intervals(l_ts)
    any(t -> any(t .∉ l_intervals) && ray.tmin < t < ray.tmax, r_ts) ||
    any(t -> any(t .∉ r_intervals) && ray.tmin < t < ray.tmax, l_ts)
end

function get_all_ts(csg::FusionCSG, ray::Ray)
    inv_ray = inv(csg.transformation) * ray
    r_ts = get_all_ts(csg.rbranch, inv_ray)
    l_ts = get_all_ts(csg.lbranch, inv_ray)
    isempty(r_ts) && return l_ts
    isempty(l_ts) && return r_ts
    r_intervals = valid_intervals(r_ts)
    l_intervals = valid_intervals(l_ts)
    r_filter = filter(t -> any(t .∉ l_intervals), r_ts)
    l_filter = filter(t -> any(t .∉ r_intervals) && t ∉ r_filter, l_ts)
    append!(r_filter, l_filter)
end


text:
"""
    pure_nash(nfg; ntofind=Inf, tol=1e-8)

Finds all pure action Nash equilibria for a normal form
game. It returns an empty array if there is no pure
action Nash.

Currently uses a brute force algorithm, but that hopefully
will change in the future.

# Arguments

- `nfg::NormalFormGame`: Instance of N-player NormalFormGame.
- `ntofind::Inf`: Maximal number of pure action Nash equilibria to be
  found; default is `prod(nfg.nums_actions)`.
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns
- `ne::Vector{NTuple{N,Int}}`: Vector of pure action Nash equilibria.
"""
function pure_nash(nfg::NormalFormGame; ntofind=prod(nfg.nums_actions),
                   tol::Real=1e-8)
    # Get number of players and their actions
    np = num_players(nfg)
    na = nfg.nums_actions

    # Holder for all NE
    ne = Array{PureActionProfile{np,Int}}(undef, 0)

    # Create counter for how many to find
    nfound = 0

    for _a in CartesianIndices(na)
        if is_nash(nfg, _a.I, tol=tol)
            push!(ne, _a.I)
            nfound = nfound + 1
        end
        nfound >= ntofind && break
    end

    return ne
end


text:
"""
    create_upper_tri(vec,V)

Create an upper triangluar matrix from a vector of the form [12, ... 1V,23,...(V-1)V]
to the form [0 12 13  ... 1V]
            [0 0  23  ... 2V]
            [...............]
            [0 0  0...(V-1)V]
# Arguments
- `vec`: vector containing values to put into the upper triangluar matrix
- `V`  : dimension of output matrix

# Returns
Upper triangluar matrix containing values of `vec`
"""
function create_upper_tri(vec::Array{T,1},V::Int64) where T <: Union{Int64,Float64}
    mat = zeros(T,V,V)
    vec2 = deepcopy(vec)
    for k = 1:V
        for l = k+1:V
            mat[k,l] = popfirst!(vec2)
        end
    end
    return mat
end

"""
    upper_triangle(matrix)

Return the upper triangle (without the diagonal) of the matrix as a vector

# Arguments
- `matrix`: matrix of which to capture the upper triangle

# Returns
Vector of upper triangluar section of `matrix`
"""
function upper_triangle(matrix::Array{T,2}) where T <: Union{Int64,Float64}
    k = 1
    ret = zeros(T,convert(Int64, round(size(matrix,1)*(size(matrix,2) - 1)/2)))
    for i in 1:size(matrix,1)
        for j in (i+1):size(matrix,2)
            ret[k] = matrix[i,j]
            k = k + 1
        end
    end
    return ret
end


text:
# one argument

@inline function eftSqr_inline{T<:StdFloat}(a::T)
    hi = a * a
    hi, fma(a, a, -hi)
end

@inline function eftIncr_inline{T<:StdFloat}(a::T)
    b = one(T)
    hi = a + b
    t = hi - a
    hi, ((a - (hi - t)) + (b - t))
end

@inline function eftDecr_inline{T<:StdFloat}(a::T)
    b = one(T)
    hi = a - b
    t = hi - a
    hi, ((a - (hi - t)) - (b + t))
end

# two arguments

@inline function eftAdd_inline{T<:StdFloat}(a::T, b::T)
  hi = a + b
  t = hi - a
  hi, ((a - (hi - t)) + (b - t))
end

@inline function eftAddGTE_inline{T<:StdFloat}(a::T, b::T)
  hi = a + b
  hi, (b - (hi - a))
end

@inline function eftSub_inline{T<:StdFloat}(a::T, b::T)
  hi = a - b
  t = hi - a
  hi, ((a - (hi - t)) - (b + t))
end

@inline function eftSubGTE_inline{T<:StdFloat}(a::T, b::T)
  hi = a - b
  hi, ((a - hi) - b)
end

@inline function eftMul_inline{T<:StdFloat}(a::T, b::T)
    hi = a * b
    hi, fma(a, b, -hi)
end

# three arguments

@inline function eftAdd_inline{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAdd_inline(b, c)
    x,u = eftAdd_inline(a, s)
    y,z = eftAdd_inline(u, t)
    x,y = eftAddGTE_inline(x, y)
    x,y,z
end

@inline function eftAddGTE_inline{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAddGTE_inline(b, c)
    x,u = eftAddGTE_inline(a, s)
    y,z = eftAddGTE_inline(u, t)
    x,y = eftAddGTE_inline(x, y)
    x,y,z
end

@inline function eftAddAs2_inline{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAdd_inline(b, c)
    x,u = eftAdd_inline(a, s)
    y   = u+t
    x,y = eftAddGTE_inline(x, y)
    x,y
end

@inline function eftMul_inline{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul_inline(a,b)
    x,y = eftMul_inline(p,c)
    z,t = eftMul_inline(e,c)
    y,z = eftAdd_inline(y,z)
    x,y = eftAdd_inline(x,y)
    eftAddGTEas3(x,y,z,t)
end

@inline function eftMulAs2_inline{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul_inline(a,b)
    x,y = eftMul_inline(p,c)
    z   = e*c
    y,z = eftAdd_inline(y,z)
    eftAdd(x,y)
end

#=
@inline function eftMulAs4_inline{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul_inline(a,b)
    x,y = eftMul_inline(p,c)
    z,t = eftMul_inline(e,c)
    y,z = eftAdd_inline(y,z)
    x,y = eftAdd_inline(x,y)
    eftAddGTE(x,y,z,t)
end
=#


text:
default(size = 0.6 .* [800, 600])
default(palette = :tol_bright)
default(dpi = 300)
default(lw = 2)
default(margin = 1mm)
FontSize = 12
default(xtickfontsize = FontSize)
default(ytickfontsize = FontSize)
default(xguidefontsize = FontSize)
default(yguidefontsize = FontSize)
default(legendfontsize = FontSize)
default(titlefontsize = FontSize)
rng = MersenneTwister(1234)

function path(θ, p = 1, ArrOut = false)
    if p == 1
        nv = size(pv1, 1)
        v = pv1
    elseif p == 2
        nv = size(pv2, 1)
        v = pv2
    else
        error("invalid path choice at path function!")
    end
    ns = nv - 1
    d = zeros(ns) #distances
    for i = 1:ns
        d[i] = EuclideanDistance(v[i, :], v[i+1, :])
    end
    td = sum(d)
    Dθ = d ./ td .* Mp
    s = 1
    thisSegBound = 0
    for i = 1:ns
        thisSegBound = sum(Dθ[1:i]) #how far this segment goes
        if θ <= thisSegBound
            s = i
            break
        end
    end
    if s != 1 #not first segment, adjust bounds and theta
        prevSegBound = sum(Dθ[1:(s-1)])
        θ -= prevSegBound
        thisSegBound -= prevSegBound
    end
    x = v[s, 1] + (v[s+1, 1] - v[s, 1]) * θ / thisSegBound
    y = v[s, 2] + (v[s+1, 2] - v[s, 2]) * θ / thisSegBound
    if θ > Mp
        x = v[end, 1]
        y = v[end, 2]
    end
    if ArrOut
        return [x, y]
    end
    return x, y
end

function EuclideanDistance(x1, x2)
    sqrt((x1[1] - x2[1])^2 + (x1[2] - x2[2])^2)
end

function uav_dynamics(x, y, vx, vy, dt, rem_power = Inf, vmax = Inf, m = 1.0)
    vx > vmax ? vmax : vx
    vy > vmax ? vmax : vy
    vx < -vmax ? -vmax : vx
    vy < -vmax ? -vmax : vy
    x = x + vx * dt
    y = y + vy * dt
    rem_power = rem_power - vx^2 * m * dt - vy^2 * m * dt - m * alpha * dt
    return x, y, rem_power
end

function simDeterministicUniform(Er = 8000)
    clearconsole()
    UASPos = [800, 450]
    LPos = [400, 550]
    ts = 0
    PNRStat = false
    N = 100
    p = 2
    dt = 1
    PrevPNR = [0.0, 0.0]
    μ = 80
    Σ = 40
    Ns = 5
    β = 0.1
    γ = 1.0
    UASPosVec = zeros(N, 2)
    anim = @animate for i = 1:N
        UASPosVec[i, :] = UASPos
        #TimeSamples = ts:120
        TimeSamples = SampleTime(μ, Σ, N, ts)
        TimeSamples = TimeSamples[:]
        PosSamples = DriverPosFunction(TimeSamples)
        elites = rankSampleMean(TimeSamples, UASPos, Ns, p)
        elite = rankSampleMean(TimeSamples, UASPos, 1, p)
        μn = sum(TimeSamples[elites]) / Ns
        μ = γ * μn + (1 - γ) * μ
        Σ =
            β * (sum((TimeSamples[elites] .- μ) .^ 2) ./ Ns + 1) + (1 - β) * Σ
        Σ = min(Σ, 100)
        OptTimeSample = TimeSamples[elite]
        v, t = DeterministicUniformMPC(
            UASPos,
            LPos,
            OptTimeSample,
            Er,
            ts,
            PrevPNR,
            PNRStat,
            p,
        )
        pp = plot()
        pp = plot!(UASPosVec[1:i, 1], UASPosVec[1:i, 2], legend = false)
        pp = drawConvexHull(TimeSamples, [1, 2], :green)
        pp = drawConvexHull(TimeSamples[elites], [1, 2], :red)
        pp = plotPlan!(UASPos, LPos, v, t, TimeSamples)
        pp = scatter!(
            path(DriverPosFunction(ts), p),
            markersize = 6.0,
            label = "",
            markercolor = :green,
        )
        pp = plot!(tickfont = Plots.font("serif", pointsize = round(12.0)))
        ts = ts + dt
        PrevPNR = v[:, 1] .* t[1] .+ UASPos
        Ed =
            m[1] * v[1, 1]^2 * t[1] +
            m[1] * v[1, 2]^2 * t[2] +
            m[2] * v[1, 3]^2 * t[3] +
            m[1] * v[2, 1]^2 * t[1] +
            m[1] * v[2, 2]^2 * t[2] +
            m[2] * v[2, 3]^2 * t[3] +
            m[1] * alpha * t[1] +
            m[1] * alpha * t[2] +
            m[2] * alpha * t[3]
        x, y, Er = uav_dynamics(
            UASPos[1],
            UASPos[2],
            v[1, 1],
            v[2, 1],
            dt,
            Er,
            vmax,
            PNRStat ? m[2] : m[1],
        )
        UASPos = [x, y]
        @show t, i, Σ, μ, Er, Ed
        if p == 2
            pv = pv2
        else
            pv = pv1
        end
        if t[1] <= 1.0 ||
           EuclideanDistance(path(DriverPosFunction(ts), p, true), pv[2, :]) <= 10
            break
        end
    end
    # proceed with rendezvous
    gif(anim, "RDV_Anim.gif", fps = 15)
end

function simIndependentDistributions(Er = 8000)
    clearconsole()
    UASPos = [800, 450]
    LPos = [400, 550]
    ts = 0
    PNRStat = false
    N = 100
    p = [1, 2]
    ptgt = 2
    dt = 1
    PrevPNR = [0.0, 0.0]
    μ = 80 * ones(2)
    Σ = 40 * ones(2)
    Ns = 5
    β = 1.0
    γ = 1.0
    UASPosVec = zeros(N, 2)
    OptTimeSample = zeros(length(p))
    anim = @animate for i = 1:N
        UASPosVec[i, :] = UASPos
        #TimeSamples = ts:120
        TimeSamples = SampleTime(μ, Σ, N, ts)
        PosSamples = DriverPosFunction(TimeSamples)
        elites = rankMultiPath(TimeSamples, UASPos, Ns, p)
        elite = rankMultiPath(TimeSamples, UASPos, 1, p)
        for j = 1:length(p)
            μn = sum(TimeSamples[elites[:, j], j]) / Ns
            μ[j] = γ * μn + (1 - γ) * μ[j]
            Σ[j] =
                β * (
                    sum((TimeSamples[elites[:, j], j] .- μ[j]) .^ 2) ./ Ns + 1
                ) + (1 - β) * Σ[j]
            Σ[j] = min(Σ[j], 100)
            OptTimeSample[j] = TimeSamples[elite[j], j]
        end
        ptgt = argmin(OptTimeSample)
        v, t = DeterministicUniformMPC(
            UASPos,
            LPos,
            OptTimeSample[2],
            Er,
            ts,
            PrevPNR,
            PNRStat,
            2,
        )
        #TODO: Fix plotting of convex hull
        pp = plot()
        pp = plot!(UASPosVec[1:i, 1], UASPosVec[1:i, 2], legend = false)
        pp = drawConvexHull(TimeSamples, [1, 2], :green)
        pp = drawConvexHull(TimeSamples[elites], [1, 2], :red)
        pp = plotPlan!(UASPos, LPos, v, t, TimeSamples)
        pp = scatter!(
            path(DriverPosFunction(ts), ptgt),
            markersize = 6.0,
            label = "",
            markercolor = :green,
        )
        pp = plot!(tickfont = Plots.font("serif", pointsize = round(12.0)))
        ts = ts + dt
        PrevPNR = v[:, 1] .* t[1] .+ UASPos
        Ed =
            m[1] * v[1, 1]^2 * t[1] +
            m[1] * v[1, 2]^2 * t[2] +
            m[2] * v[1, 3]^2 * t[3] +
            m[1] * v[2, 1]^2 * t[1] +
            m[1] * v[2, 2]^2 * t[2] +
            m[2] * v[2, 3]^2 * t[3] +
            m[1] * alpha * t[1] +
            m[1] * alpha * t[2] +
            m[2] * alpha * t[3]
        x, y, Er = uav_dynamics(
            UASPos[1],
            UASPos[2],
            v[1, 1],
            v[2, 1],
            dt,
            Er,
            vmax,
            PNRStat ? m[2] : m[1],
        )
        UASPos = [x, y]
        @show t, i, Σ, μ, Er, Ed
        if ptgt == 2
            pv = pv2
        else
            pv = pv1
        end
        if t[1] <= 1.0 ||
           EuclideanDistance(
            path(DriverPosFunction(ts), ptgt, true),
            pv[2, :],
        ) <= 10
            break
        end
    end
    # proceed with rendezvous
    gif(anim, "RDV_Anim.gif", fps = 15)
end

max_stars_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_issues_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_forks_repo_name: gbarsih/Multi-Path-Safe-Rendezvous


text:
# module SQLite_Octo.SQL

function query(q::QueryBuilders.Queryable)
    SQLite.query(q.adapter.db, repr(q))
end


text:
## Shapefiles functions useful for mapping

# Download worldmap shapefile
function download_shapefile(res)
    @assert res ∈ [50, 100]
    dir =
        "https://github.com/nvkelso/natural-earth-vector/" * "raw/master/$(res)m_physical/"

    fn = "ne_$(res)m_land.shp"
    return run(`wget $dir/$fn -P ./assets/`)
end

# Read worlmap shapefile
function worldshape(res)
    handle = open("./assets/ne_$(res)m_land.shp", "r") do io
        read(io, Shapefile.Handle)
    end
    return handle
end

# Clip worldmap shapefile
function clip(s::Shapefile.Handle, l::SimpleSDMLayer)
    return filter(x -> isin(x, l), s.shapes)
end

# Check layer bounds are outside polygon
function isin(p::Shapefile.Polygon, l::SimpleSDMLayer)
    out = false
    for xy in p.points
        xy.x < l.right && (out = true)
        xy.x > l.left && (out = true)
        xy.y < l.top && (out = true)
        xy.y > l.bottom && (out = true)
    end
    return out
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
function equalvalued(X::NTuple)
    n = length(X)
    allequal = true
    i = 2
    while allequal && i <= n
        allequal = X[i] == X[i-1]
        i += 1
    end #while
    return allequal
end #equalvalued

sizes(As::AxisArray...) = tuple(zip(map(a -> map(length, Base.axes(a)), As)...)...)
matchingdims(As::Tuple{Vararg{AxisArray}}) = all(equalvalued, sizes(As...))
matchingdimsexcept(As::Tuple{Vararg{AxisArray}}, n::Int) = all(equalvalued, sizes(As...)[[1:n-1; n+1:end]])

Base.cat(A1::AxisArray{T}, As::AxisArray{T}...; dims) where {T} = _cat(dims, A1, As...)
_cat(::Val{n}, As...) where {n} = _cat(n, As...)

@inline function _cat(n::Integer, As...)
    if n <= ndims(As[1])
        matchingdimsexcept(As, n) || error("All non-concatenated axes must be identically-valued")
        newaxis = Axis{axisnames(As[1])[n]}(vcat(map(A -> A.axes[n].val, As)...))
        checkaxis(newaxis)
        return AxisArray(cat(map(A->A.data, As)..., dims=n), (As[1].axes[1:n-1]..., newaxis, As[1].axes[n+1:end]...))
    else
        matchingdims(As) || error("All axes must be identically-valued")
        return AxisArray(cat(map(A->A.data, As)..., dims=n), As[1].axes)
    end #if
end

function axismerge(method::Symbol, axes::Axis{name,T}...) where {name,T}

    axisvals = if method == :inner
        intersect(axisvalues(axes...)...)
    elseif method == :left
        axisvalues(axes[1])[1]
    elseif method == :right
        axisvalues(axes[end])[1]
    elseif method == :outer
        union(axisvalues(axes...)...)
    else
        error("Join method must be one of :inner, :left, :right, :outer")
    end #if

    isa(axistrait(axisvals), Dimensional) && sort!(axisvals)

    return Axis{name}(collect(axisvals))

end

function indexmappings(oldaxes::NTuple{N,Axis}, newaxes::NTuple{N,Axis}) where N
    oldvals = axisvalues(oldaxes...)
    newvals = axisvalues(newaxes...)
    return collect(zip(indexmapping.(oldvals, newvals)...))
end

function indexmapping(old::AbstractVector, new::AbstractVector)

    before = Int[]
    after = Int[]

    oldperm = sortperm(old)
    newperm = sortperm(new)

    oldsorted = old[oldperm]
    newsorted = new[newperm]

    oldlength = length(old)
    newlength = length(new)

    oi = ni = 1

    while oi <= oldlength && ni <= newlength

        oldval = oldsorted[oi]
        newval = newsorted[ni]

        if oldval == newval
            push!(before, oldperm[oi])
            push!(after, newperm[ni])
            oi += 1
            ni += 1
        elseif oldval < newval
            oi += 1
        else
            ni += 1
        end

    end

    return before, after

end

"""
    merge(As::AxisArray...)

Combines AxisArrays with matching axis names into a single AxisArray spanning all of the axis values of the inputs. If a coordinate is defined in more than ones of the inputs, it takes its value from last input in which it appears. If a coordinate in the output array is not defined in any of the input arrays, it takes the value of the optional `fillvalue` keyword argument (default zero).
"""
function Base.merge(As::AxisArray{T,N,D,Ax}...; fillvalue::T=zero(T)) where {T,N,D,Ax}

    resultaxes = map(as -> axismerge(:outer, as...), map(tuple, axes.(As)...))
    resultdata = fill(fillvalue, length.(resultaxes)...)
    result = AxisArray(resultdata, resultaxes...)

    for A in As
        before_idxs, after_idxs = indexmappings(A.axes, result.axes)
        result.data[after_idxs...] = A.data[before_idxs...]
    end

    return result

end #merge

"""
    join(As::AxisArray...)

Combines AxisArrays with matching axis names into a single AxisArray. Unlike `merge`, the inputs are joined along a newly created axis (optionally specified with the `newaxis` keyword argument).  The `method` keyword argument can be used to specify the join type:

`:inner` - keep only those array values at axis values common to all AxisArrays to be joined
`:left` - keep only those array values at axis values present in the first AxisArray passed
`:right` - keep only those array values at axis values present in the last AxisArray passed
`:outer` (default) - keep all array values: create an AxisArray spanning all of the input axis values

If an array value in the output array is not defined in any of the input arrays (i.e. in the case of a left, right, or outer join), it takes the value of the optional `fillvalue` keyword argument (default zero).
"""
function Base.join(As::AxisArray{T,N,D,Ax}...; fillvalue::T=zero(T),
                   newaxis::Axis=_default_axis(1:length(As), ndims(As[1])+1),
                   method::Symbol=:outer) where {T,N,D,Ax}

    prejoin_resultaxes = map(as -> axismerge(method, as...), map(tuple, axes.(As)...))

    resultaxes = (prejoin_resultaxes..., newaxis)
    resultdata = fill(fillvalue, length.(resultaxes)...)
    result = AxisArray(resultdata, resultaxes...)

    for (i, A) in enumerate(As)
        before_idxs, after_idxs = indexmappings(A.axes, prejoin_resultaxes)
        result.data[(after_idxs..., i)...] = A.data[before_idxs...]
    end #for

    return result

end #join

function _collapse_array_axes(array_name, array_axes...)
    ((array_name, (idx isa Tuple ? idx : (idx,))...) for idx in Iterators.product((Ax.val for Ax in array_axes)...))
end

function _collapse_axes(array_names, array_axes)
    collect(Iterators.flatten(map(array_names, array_axes) do tup_name, tup_array_axes
        _collapse_array_axes(tup_name, tup_array_axes...)
    end))
end

function _splitall(::Val{N}, As...) where N
    tuple((Base.IteratorsMD.split(A, Val(N)) for A in As)...)
end

function _reshapeall(::Val{N}, As...) where N
    tuple((reshape(A, Val(N)) for A in As)...)
end

function _check_common_axes(common_axis_tuple)
    if !all(axisname(first(common_axis_tuple)) .=== axisname.(common_axis_tuple[2:end]))
        throw(ArgumentError("Leading common axes must have the same name in each array"))
    end

    return nothing
end

function _collapsed_axis_eltype(LType, trailing_axes)
    eltypes = map(trailing_axes) do array_trailing_axes
        Tuple{LType, eltype.(array_trailing_axes)...}
    end

    return typejoin(eltypes...)
end

function collapse(::Val{N}, As::Vararg{AxisArray, AN}) where {N, AN}
    collapse(Val(N), ntuple(identity, AN), As...)
end

function collapse(::Val{N}, ::Type{NewArrayType}, As::Vararg{AxisArray, AN}) where {N, AN, NewArrayType<:AbstractArray}
    collapse(Val(N), NewArrayType, ntuple(identity, AN), As...)
end

@generated function collapse(::Val{N}, labels::NTuple{AN, LType}, As::Vararg{AxisArray, AN}) where {N, AN, LType}
    collapsed_dim_int = Int(N) + 1
    new_eltype = Base.promote_eltype(As...)

    quote
        collapse(Val(N), Array{$new_eltype, $collapsed_dim_int}, labels, As...)
    end
end

"""
    collapse(::Val{N}, As::AxisArray...) -> AxisArray
    collapse(::Val{N}, labels::Tuple, As::AxisArray...) -> AxisArray
    collapse(::Val{N}, ::Type{NewArrayType}, As::AxisArray...) -> AxisArray
    collapse(::Val{N}, ::Type{NewArrayType}, labels::Tuple, As::AxisArray...) -> AxisArray

Collapses `AxisArray`s with `N` equal leading axes into a single `AxisArray`.
All additional axes in any of the arrays are collapsed into a single additional
axis of type `Axis{:collapsed, CategoricalVector{Tuple}}`.

### Arguments

* `::Val{N}`: the greatest common dimension to share between all input
                    arrays. The remaining axes are collapsed. All `N` axes must be common
                    to each input array, at the same dimension. Values from `0` up to the
                    minimum number of dimensions across all input arrays are allowed.
* `labels::Tuple`: (optional) an index for each array in `As` used as the leading element in
                   the index tuples in the `:collapsed` axis. Defaults to `1:length(As)`.
* `::Type{NewArrayType<:AbstractArray{_, N+1}}`: (optional) the desired underlying array
                                                 type for the returned `AxisArray`.
* `As::AxisArray...`: `AxisArray`s to be collapsed together.

### Examples

```
julia> price_data = AxisArray(rand(10), Axis{:time}(Date(2016,01,01):Day(1):Date(2016,01,10)))
1-dimensional AxisArray{Float64,1,...} with axes:
    :time, 2016-01-01:1 day:2016-01-10
And data, a 10-element Array{Float64,1}:
 0.885014
 0.418562
 0.609344
 0.72221
 0.43656
 0.840304
 0.455337
 0.65954
 0.393801
 0.260207

julia> size_data = AxisArray(rand(10,2), Axis{:time}(Date(2016,01,01):Day(1):Date(2016,01,10)), Axis{:measure}([:area, :volume]))
2-dimensional AxisArray{Float64,2,...} with axes:
    :time, 2016-01-01:1 day:2016-01-10
    :measure, Symbol[:area, :volume]
And data, a 10×2 Array{Float64,2}:
 0.159434     0.456992
 0.344521     0.374623
 0.522077     0.313256
 0.994697     0.320953
 0.95104      0.900526
 0.921854     0.729311
 0.000922581  0.148822
 0.449128     0.761714
 0.650277     0.135061
 0.688773     0.513845

julia> collapsed = collapse(Val(1), (:price, :size), price_data, size_data)
2-dimensional AxisArray{Float64,2,...} with axes:
    :time, 2016-01-01:1 day:2016-01-10
    :collapsed, Tuple{Symbol,Vararg{Symbol,N} where N}[(:price,), (:size, :area), (:size, :volume)]
And data, a 10×3 Array{Float64,2}:
 0.885014  0.159434     0.456992
 0.418562  0.344521     0.374623
 0.609344  0.522077     0.313256
 0.72221   0.994697     0.320953
 0.43656   0.95104      0.900526
 0.840304  0.921854     0.729311
 0.455337  0.000922581  0.148822
 0.65954   0.449128     0.761714
 0.393801  0.650277     0.135061
 0.260207  0.688773     0.513845

julia> collapsed[Axis{:collapsed}(:size)] == size_data
true
```

"""
@generated function collapse(::Val{N},
                             ::Type{NewArrayType},
                             labels::NTuple{AN, LType},
                             As::Vararg{AxisArray, AN}) where {N, AN, LType, NewArrayType<:AbstractArray}
    if N < 0
        throw(ArgumentError("collapse dimension N must be at least 0"))
    end

    if N > minimum(ndims.(As))
        throw(ArgumentError(
            """
            collapse dimension N must not be greater than the maximum number of dimensions
            across all input arrays
            """
        ))
    end

    collapsed_dim = Val(N + 1)
    collapsed_dim_int = Int(N) + 1

    common_axes, trailing_axes = zip(_splitall(Val(N), axisparams.(As)...)...)

    foreach(_check_common_axes, zip(common_axes...))

    new_common_axes = first(common_axes)
    collapsed_axis_eltype = _collapsed_axis_eltype(LType, trailing_axes)
    collapsed_axis_type = CategoricalVector{collapsed_axis_eltype, Vector{collapsed_axis_eltype}}

    new_axes_type = Tuple{new_common_axes..., Axis{:collapsed, collapsed_axis_type}}
    new_eltype = Base.promote_eltype(As...)

    quote
        common_axes, trailing_axes = zip(_splitall(Val(N), axes.(As)...)...)

        for common_axis_tuple in zip(common_axes...)
            if !isempty(common_axis_tuple)
                for common_axis in common_axis_tuple[2:end]
                    if !all(axisvalues(common_axis) .== axisvalues(common_axis_tuple[1]))
                        throw(ArgumentError(
                            """
                            Leading common axes must be identical across
                            all input arrays"""
                        ))
                    end
                end
            end
        end

        array_data = cat(_reshapeall($collapsed_dim, As...)..., dims=$collapsed_dim)

        axis_array_type = AxisArray{
            $new_eltype,
            $collapsed_dim_int,
            $NewArrayType,
            $new_axes_type
        }

        new_axes = (
            first(common_axes)...,
            Axis{:collapsed, $collapsed_axis_type}($collapsed_axis_type(_collapse_axes(labels, trailing_axes))),
        )

        return axis_array_type(array_data, new_axes)
    end
end


text: include("testGetData.jl")


text:
using Test
using Dates

include("../src/PortableGameNotation.jl")
using .PortableGameNotation

g = readpgn("test.pgn")

@testset "1" begin
    @test date(g[1]) == Date(2015,7,3)
    @test site(g[1]) == "Danzhou CHN"
    @test event(g[1]) == "6th Hainan Danzhou"
    @test result(g[1]) == "1-0"
    @test intresult(g[1]) == 1
    @test white(g[1]) == "Wei Yi"
    @test black(g[1]) == "Lazaro Bruzon Batista"
    @test whiteelo(g[1]) == 2724
    @test blackelo(g[1]) == 2669
    @test eco(g[1]) == "B40"
    @test eventdate(g[1]) == Date(2015,7,2)
    @test plycount(g[1]) == 71
    #@test movetext(g[1]) == """1. e4 c5 2. Nf3 e6 3. Nc3 a6 4. Be2 Nc6 5. d4 cxd4 6. Nxd4 Qc7 7. O-O Nf6 8. Be3 Be7 9. f4 d6 10. Kh1 O-O 11. Qe1 Nxd4 12. Bxd4 b5 13. Qg3 Bb7 14. a3 Rad8 15. Rae1 Rd7 16. Bd3 Qd8 17. Qh3 g6 18. f5 e5 19. Be3 Re8 20. fxg6 hxg6 21. Nd5 Nxd5 22. Rxf7 Kxf7 23. Qh7+ Ke6 24. exd5+ Kxd5 25. Be4+ Kxe4 26. Qf7 Bf6 27. Bd2+ Kd4 28. Be3+ Ke4 29. Qb3 Kf5 30. Rf1+ Kg4 31. Qd3 Bxg2+ 32. Kxg2 Qa8+ 33. Kg1 Bg5 34. Qe2+ Kh4 35. Bf2+ Kh3 36. Be1 1-0"""
    @test length(g[1]) == 36
    @test whiteev(g[1]) == 0.5784967523447427
    @test blackev(g[1]) == 0.42150324765525726
    @test whitescore(g[1]) == 1
    @test blackscore(g[1]) == 0
    @test whiteperfelo(g[1]) == 3069
    @test blackperfelo(g[1]) == 2324
    @test isdecisive(g[1]) == true
    @test PortableGameNotation.validate(g[1])
end

@testset "2" begin
    @test date(g[2]) == Date(2016,4,29)
    @test site(g[2]) == "St. Louis, MO USA"
    @test event(g[2]) == "Ultimate Blitz Challenge"
    @test result(g[2]) == "1-0"
    @test intresult(g[2]) == 1
    @test white(g[2]) == "Wesley So"
    @test black(g[2]) == "Garry Kasparov"
    @test whiteelo(g[2]) == 2773
    @test blackelo(g[2]) == 2812
    @test eventdate(g[2]) == Date(2016,4,28)
    @test eco(g[2]) == "A41"
    @test plycount(g[2]) == 49
    #@test movetext(g[2]) == """1. Nf3 g6 2. e4 Bg7 3. d4 d6 4. c4 Bg4 5. Be2 Nc6 6. Nbd2 e5 7. d5 Nce7 8. h3 Bd7 9. c5 dxc5 10. Nc4 f6 11. d6 Nc8 12. Be3 b6 13. O-O Bc6 14. dxc7 Qxc7 15. b4 cxb4 16. Rc1 Nge7 17. Qb3 h6 18. Rfd1 b5 19. Ncxe5 fxe5 20. Bxb5 Rb8 21. Ba4 Qb7 22. Rxc6 Nxc6 23. Qe6+ Ne7 24. Bc5 Rc8 25. Bxe7 1-0"""
    @test length(g[2]) == 25
    @test whiteev(g[2]) == 0.4441090388831469
    @test blackev(g[2]) == 0.5558909611168531
    @test whitescore(g[2]) == 1
    @test blackscore(g[2]) == 0
    @test whiteperfelo(g[2]) == 3212
    @test blackperfelo(g[2]) == 2373
    @test isdecisive(g[2]) == true
    @test PortableGameNotation.validate(g[2])
end

games = readpgn("testfen.pgn")

for g in games
    println(g)
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/variable.jl
# Testing for Variable
#############################################################################
using JuMP, FactCheck

facts("[variable] constructors") do
    # Constructors
    mcon = Model()
    @variable(mcon, nobounds)
    @variable(mcon, lbonly >= 0)
    @variable(mcon, ubonly <= 1)
    @variable(mcon, 0 <= bothb <= 1)
    @variable(mcon, 0 <= onerange[-5:5] <= 10)
    @variable(mcon, onerangeub[-7:1] <= 10, Int)
    @variable(mcon, manyrangelb[0:1,10:20,1:1] >= 2)
    @fact getlowerbound(manyrangelb[0,15,1]) --> 2
    s = ["Green","Blue"]
    @variable(mcon, x[i=-10:10,s] <= 5.5, Int, start=i+1)
    @fact getupperbound(x[-4,"Green"]) --> 5.5
    @fact getvalue(x[-3,"Blue"]) --> -2
    @fact isequal(getvariable(mcon, :lbonly),lbonly) --> true
    @fact isequal(getvariable(mcon, :ubonly),ubonly) --> true
    @fact isequal(getvariable(mcon, :onerangeub)[-7],onerangeub[-7]) --> true
    @variable(mcon, lbonly)
    @fact_throws ErrorException getvariable(mcon, :lbonly)
    @fact_throws ErrorException getvariable(mcon, :foo)
    d = Dict()
    @variable(mcon, d["bar"][1:10] == 1)
    @fact getvalue(d["bar"][1]) --> 1
    @fact typeof(zero(nobounds)) --> AffExpr
    @fact typeof(one(nobounds)) --> AffExpr
end

facts("[variable] get and set bounds") do
    m = Model()
    @variable(m, 0 <= x <= 2)
    @fact getlowerbound(x) --> 0
    @fact getupperbound(x) --> 2
    setlowerbound(x, 1)
    @fact getlowerbound(x) --> 1
    setupperbound(x, 3)
    @fact getupperbound(x) --> 3
    @variable(m, y, Bin)
    @fact getlowerbound(y) --> 0
    @fact getupperbound(y) --> 1
    @variable(m, 0 <= y <= 1, Bin)
    @fact getlowerbound(y) --> 0
    @fact getupperbound(y) --> 1
    @variable(m, fixedvar == 2)
    @fact getvalue(fixedvar) --> 2
    @fact getlowerbound(fixedvar) --> 2
    @fact getupperbound(fixedvar) --> 2
    setvalue(fixedvar, 5)
    @fact getvalue(fixedvar) --> 5
    @fact getlowerbound(fixedvar) --> 5
    @fact getupperbound(fixedvar) --> 5
end

facts("[variable] get and set values") do
    m = Model()
    @variable(m, x[1:3])
    x0 = collect(1:3)
    setvalue(x, x0)
    @fact getvalue(x) --> x0
    @fact getvalue([x[1],x[2],x[3]]) --> x0

    @variable(m, y[1:3,1:2])
    @fact_throws DimensionMismatch setvalue(y, collect(1:6))
end

facts("[variable] get and set category") do
    m = Model()
    @variable(m, x[1:3])
    setcategory(x[2], :Int)
    @fact getcategory(x[3]) --> :Cont
    @fact getcategory(x[2]) --> :Int
end

facts("[variable] repeated elements in index set (issue #199)") do
    repeatmod = Model()
    s = [:x,:x,:y]
    @variable(repeatmod, x[s])
    @fact MathProgBase.numvar(repeatmod) --> 3
end

facts("[variable] condition in indexing") do
    fa = repl[:for_all]
    inset, dots = repl[:in], repl[:dots]
    condmod = Model()
    @variable(condmod, x[i=1:10; iseven(i)])
    @variable(condmod, y[j=1:10,k=3:2:9; isodd(j+k) && k <= 8])
    @fact length(x.tupledict) --> 5
    @fact length(y.tupledict) --> 15
    @fact string(condmod) --> "Min 0\nSubject to\n x[i] free $fa i $inset {1,2,$dots,9,10} s.t. iseven(i)\n y[j,k] free $fa j $inset {1,2,$dots,9,10}, k $inset {3,5,7,9} s.t. isodd(j + k) and k <= 8\n"
end

facts("[variable] @variable returning Array{Variable}") do
    m = Model()
    @variable(m, x[1:3,1:4,1:2])
    @variable(m, y[1:0])
    @variable(m, z[1:4])

    @fact typeof(x) --> Array{Variable,3}
    @fact typeof(y) --> Array{Variable,1}
    @fact typeof(z) --> Array{Variable,1}

    @fact typeof(getvalue(x)) --> Array{Float64,3}
    @fact typeof(getvalue(y)) --> Array{Float64,1}
    @fact typeof(getvalue(z)) --> Array{Float64,1}
end

facts("[variable] getvalue on empty things") do
    m = Model()
    @variable(m, x[1:4,  1:0,1:3])   # Array{Variable}
    @variable(m, y[1:4,  2:1,1:3]) # JuMPArray
    @variable(m, z[1:4,Set(),1:3]) # JuMPDict

    @fact getvalue(x) --> Array(Float64, 4, 0, 3)
    @fact typeof(getvalue(y)) <: JuMP.JuMPArray{Float64} --> true
    @fact JuMP.size(getvalue(y)) --> (4,0,3)
    @fact typeof(getvalue(z)) --> JuMP.JuMPArray{Float64,3,Tuple{UnitRange{Int},Set{Any},UnitRange{Int}}}
    @fact length(getvalue(z)) --> 0
end

# Slices three-dimensional JuMPContainer x[I,J,K]
# I,J,K can be singletons, ranges, colons, etc.
function sliceof(x, I, J, K)
    y = Array(Variable, length(I), length(J), length(K))

    ii = 1
    jj = 1
    kk = 1
    for i in I
        for j in J
            for k in K
                y[ii,jj,kk] = x[i,j,k]
                kk += 1
            end
            jj += 1
            kk = 1
        end
        ii += 1
        jj = 1
    end
    if VERSION >= v"0.5.0-dev+1195"
        idx = [length(I)==1, length(J)==1, length(K)==1]
    else
        idx = [false, false, length(K)==1]
        if idx[3] == true && length(J) == 1
            idx[2] = true
        end
    end
    squeeze(y, tuple(find(idx)...))
end

facts("[variable] Slices of JuMPArray (#684)") do
    m = Model()
    @variable(m, x[1:3, 1:4,1:2])
    @variable(m, y[1:3,-1:2,3:4])
    @variable(m, z[1:3,-1:2:4,3:4])
    @variable(m, w[1:3,-1:2,[:red,"blue"]])

    @fact x[:] --> vec(sliceof(x, 1:3, 1:4, 1:2))
    @fact x[:,:,:] --> sliceof(x, 1:3, 1:4, 1:2)
    @fact x[1,:,:] --> sliceof(x, 1, 1:4, 1:2)
    @fact x[1,:,2] --> sliceof(x, 1, 1:4, 2)
    @fact_throws x[1,:,3]
    @fact x[1:2,:,:] --> sliceof(x, 1:2, 1:4, 1:2)
    @fact x[1:2,:,2] --> sliceof(x, 1:2, 1:4, 2)
    @fact x[1:2,:,1:2] --> sliceof(x, 1:2, 1:4, 1:2)
    @fact_throws x[1:2,:,1:3]

    @fact y[:] --> vec(sliceof(y, 1:3, -1:2, 3:4))
    @fact y[:,:,:] --> sliceof(y, 1:3, -1:2, 3:4)
    @fact y[1,:,:] --> sliceof(y, 1, -1:2, 3:4)
    @fact y[1,:,4] --> sliceof(y, 1, -1:2, 4)
    @fact_throws y[1,:,5]
    @fact y[1:2,:,:] --> sliceof(y, 1:2, -1:2, 3:4)
    @fact y[1:2,:,4] --> sliceof(y, 1:2, -1:2, 4)
    @fact y[1:2,:,3:4] --> sliceof(y, 1:2, -1:2, 3:4)
    @fact_throws y[1:2,:,1:3]

    @fact z[:] --> vec(sliceof(z, 1:3, -1:2:4, 3:4))
    @fact z[:,1,:] --> sliceof(z, 1:3, 1, 3:4)
    @fact z[1,1,:] --> sliceof(z, 1, 1, 3:4)
    @fact_throws z[:,5,3]
    @fact z[1:2,1,:] --> sliceof(z, 1:2, 1, 3:4)
    @fact z[1:2,1,4] --> sliceof(z, 1:2, 1, 4)
    @fact z[1:2,1,3:4] --> sliceof(z, 1:2, 1, 3:4)
    @fact_throws z[1:2,1,1:3]

    @fact w[:] --> vec(sliceof(w, 1:3, -1:2, [:red,"blue"]))
    @fact_throws w[:,:,:]
    @fact w[1,:,"blue"] --> sliceof(w, 1, -1:2, ["blue"])
    @fact w[1,:,:red] --> sliceof(w, 1, -1:2, [:red])
    @fact_throws w[1,:,"green"]
    @fact w[1:2,:,"blue"] --> sliceof(w, 1:2, -1:2, ["blue"])
    @fact_throws w[1:2,:,[:red,"blue"]]
end

facts("[variable] Can't use end for indexing a JuMPContainer") do
    m = Model()
    @variable(m, x[0:2,1:4])
    @variable(m, y[i=1:4,j=1:4;true])
    @variable(m, z[0:2])
    @fact_throws x[end,1]
    @fact_throws x[end-1]
    @fact_throws x[0,end-1]
    @fact_throws y[end,end-1]
    @fact_throws y[end,1]
    @fact_throws z[end]
end


text:
import MathOptInterface
import MathOptInterface: Utilities

const MOI = MathOptInterface
const MOIU = MOI.Utilities

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# function aliases
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}

# set aliases
const Bounds{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T}
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::BaronModel
    nlp_block_data::Union{Nothing, MOI.NLPBlockData}
end

Optimizer(; options...) = Optimizer(BaronModel(; options...), nothing)

MOI.get(model::Optimizer, ::MOI.SolverName) = "BARON"

# empty
function MOI.is_empty(model::Optimizer)
    BARON.is_empty(model.inner) && model.nlp_block_data === nothing
end

function MOI.empty!(model::Optimizer)
    model.inner = BaronModel(; ((Symbol(key), val) for (key, val) in model.inner.options)...)
    model.nlp_block_data = nothing
    return
end

# copy
# MOIU.supports_default_copy_to(::Optimizer, copy_names::Bool) = !copy_names
MOI.supports_incremental_interface(::Optimizer) = true
function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.default_copy_to(dest, src; kws...)
end

# optimize
function MOI.optimize!(model::Optimizer)
    if !IS_SOLVER_SET
        error(("""BARON.jl was not built correctly.
                 Set the environment variable `BARON_EXEC` and run `using Pkg; Pkg.build("BARON")`."""))
    end
    write_bar_file(model.inner)
    if model.inner.print_input_file
        println("\nBARON input file: $(model.inner.problem_file_name)\n")
        println(read(model.inner.problem_file_name, String))
    end
    try
        run(`$baron_exec $(model.inner.problem_file_name)`)
    catch e
        println("$e")
        println(read(model.inner.problem_file_name, String))
        error("failed to call BARON exec $baron_exec")
    end
    read_results(model.inner)
end

# RawOptimizerAttribute
MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true
function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    model.inner.options[param.name] = value
    return
end
function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    return get(model.inner.options, param.name) do
        throw(ErrorException("Requested parameter $(param.name) is not set."))
    end
end

# TimeLimitSec
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, val::Real)
    model.inner.options["MaxTime"] = Float64(val)
    return
end

# BARON's default time limit is 1000 seconds.
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.inner.options, "MaxTime", 1000.0)
end

struct PrintInputFile <: MOI.AbstractOptimizerAttribute end
function MOI.set(model::Optimizer, ::PrintInputFile, val::Bool)
    model.inner.print_input_file = val
    return
end
function MOI.get(model::Optimizer, ::PrintInputFile)
    model.inner.print_input_file
end

include(joinpath("moi", "util.jl"))
include(joinpath("moi", "variables.jl"))
include(joinpath("moi", "constraints.jl"))
include(joinpath("moi", "objective.jl"))
include(joinpath("moi", "nlp.jl"))
include(joinpath("moi", "results.jl"))


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

module BuildAGF

using Pkg.Artifacts

# paths for AGFFileReader source file builds
const GLASSCAT_DIR = @__DIR__ # contains AGFFileReader.jl (pre-existing)
const SOURCES_PATH = joinpath(GLASSCAT_DIR, "sources.txt")

include("generate.jl")

include(joinpath("..", "src", "constants.jl"))
# include(joinpath(GLASSCAT_DIR, "GlassTypes.jl"))
include(joinpath(GLASSCAT_DIR, "sources.jl"))
include(joinpath(GLASSCAT_DIR, "generate.jl"))

artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
agf_hash = artifact_hash("AGF", artifact_toml)

# If the name was not bound, or the hash it was bound to does not exist, create it!
if agf_hash === nothing || !artifact_exists(agf_hash)
    # create_artifact() returns the content-hash of the artifact directory once we're finished creating it

    agf_hash = create_artifact() do artifact_dir
        agf_dir = joinpath(artifact_dir, "agf") # contains SCHOTT.agf, SUMITA.agf, etc.
        mkdir(agf_dir)
        jl_dir = joinpath(artifact_dir, "jl") # contains AGFGlasscat.jl, SCHOTT.jl, etc.
        mkdir(jl_dir)
       #run agf generation code to create a glass catalog directory in tempdir
       #then copy this to the artifact location
       # Build/verify a source directory using information from sources.txt
        sources = split.(readlines(SOURCES_PATH))
        verify_sources!(sources, agf_dir)
        verified_source_names = first.(sources)

        # Use verified sources to generate required .jl files
        @info "Using sources: $(join(verified_source_names, ", ", " and "))"
        @info "Creating glass catalog artifacts in directory $artifact_dir with subdirectories $jl_dir $agf_dir"
        generate_jls(verified_source_names, artifact_dir, jl_dir, agf_dir)
        cp(joinpath(@__DIR__,GlassTypes.jl),artifact_dir)
    end

    # Now bind that hash within our `Artifacts.toml`.  `force = true` means that if it already exists,
    # just overwrite with the new content-hash.  Unless the source files change, we do not expect
    # the content hash to change, so this should not cause unnecessary version control churn.
    bind_artifact!(artifact_toml, "AGF", agf_hash,force=true)
end

end #module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
module DDL
using HTTP, ZipFile, Tar

end


text:
export collectLatest

import Base: show

"""
    collectLatest(sources::S, mappingFn::F = copy) where { S, F }
    collectLatest(::Type{T}, ::Type{R}, sources::S, mappingFn::F = copy)

Collects values from multible Observables and emits it in one single array every time each inner Observable has a new value.
Reemits errors from inner observables. Completes when all inner observables completes.

# Arguments
- `sources`: input sources
- `mappingFn`: optional mappingFn applied to an array of emited values, `copy` by default, should return a Vector

# Optional arguments
- `::Type{T}`: optional type of emmiting values of inner observables
- `::Type{R}`: optional return type after applying `mappingFn` to a vector of values

# Examples
```jldoctest
using Rocket

collected = collectLatest([ of(1), from([ 1, 2 ]) ])

subscribe!(collected, logger())
;

# output

[LogActor] Data: [1, 1]
[LogActor] Data: [1, 2]
[LogActor] Completed
```

See also: [`Subscribable`](@ref), [`subscribe!`](@ref), [`combineLatest`](@ref)
"""
function collectLatest(sources::S, mappingFn::F = copy) where { S, F }
    T = union_type(sources)
    R = similar_typeof(sources, T)
    return CollectLatestObservable{T, S, R, F}(sources, mappingFn)
end

collectLatest(::Type{T}, ::Type{R}, sources::S, mappingFn::F = copy) where { T, R, S, F } = CollectLatestObservable{T, S, R, F}(sources, mappingFn)

##

struct CollectLatestObservableWrapper{L, A, S, B, T, F}
    actor   :: A
    storage :: S

    cstatus       :: B # Completion status
    vstatus       :: B # Values status
    ustatus       :: B # Updates status
    subscriptions :: T
    mappingFn     :: F

    CollectLatestObservableWrapper{L, A, S, B, T, F}(actor::A, storage::S, cstatus::B, vstatus::B, ustatus::B, subscriptions::T, mappingFn::F) where {L, A, S, B, T, F} = begin
        return new(actor, storage, cstatus, vstatus, ustatus, subscriptions, mappingFn)
    end
end

function CollectLatestObservableWrapper(::Type{L}, actor::A, storage::S, mappingFn::F) where { L, A, S, F }
    nsize         = size(storage)
    cstatus       = falses(nsize)
    vstatus       = falses(nsize)
    ustatus       = falses(nsize)
    subscriptions = fill!(similar(storage, Teardown), voidTeardown)
    return CollectLatestObservableWrapper{L, A, S, typeof(cstatus), typeof(subscriptions), F}(actor, storage, cstatus, vstatus, ustatus, subscriptions, mappingFn)
end

cstatus(wrapper::CollectLatestObservableWrapper, index::CartesianIndex) = @inbounds wrapper.cstatus[index]
vstatus(wrapper::CollectLatestObservableWrapper, index::CartesianIndex) = @inbounds wrapper.vstatus[index]
ustatus(wrapper::CollectLatestObservableWrapper, index::CartesianIndex) = @inbounds wrapper.ustatus[index]

dispose(wrapper::CollectLatestObservableWrapper) = begin fill!(wrapper.cstatus, true); foreach(s -> unsubscribe!(s), wrapper.subscriptions) end

struct CollectLatestObservableInnerActor{L, I <: CartesianIndex, W} <: Actor{L}
    index   :: I
    wrapper :: W
end

Base.show(io::IO, ::CollectLatestObservableInnerActor{L}) where L = print(io, "CollectedObservableInnerActor($L)")

on_next!(actor::CollectLatestObservableInnerActor{L}, data::L) where L = next_received!(actor.wrapper, data, actor.index)
on_error!(actor::CollectLatestObservableInnerActor, err)               = error_received!(actor.wrapper, err, actor.index)
on_complete!(actor::CollectLatestObservableInnerActor)                 = complete_received!(actor.wrapper, actor.index)

function next_received!(wrapper::CollectLatestObservableWrapper, data, index::CartesianIndex)
    @inbounds wrapper.storage[index] = data
    @inbounds wrapper.vstatus[index] = true
    @inbounds wrapper.ustatus[index] = true
    if all(wrapper.vstatus) && !all(wrapper.cstatus)
        unsafe_copyto!(wrapper.vstatus, 1, wrapper.cstatus, 1, length(wrapper.vstatus))
        next!(wrapper.actor, wrapper.mappingFn(wrapper.storage))
    end
end

function error_received!(wrapper::CollectLatestObservableWrapper, err, index::CartesianIndex)
    if !(@inbounds wrapper.cstatus[index])
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::CollectLatestObservableWrapper, index::CartesianIndex)
    if !all(wrapper.cstatus)
        @inbounds wrapper.cstatus[index] = true
        if ustatus(wrapper, index)
            @inbounds wrapper.vstatus[index] = true
        end
        if all(wrapper.cstatus) || (@inbounds wrapper.vstatus[index] === false)
            dispose(wrapper)
            complete!(wrapper.actor)
        end
    end
end

##

@subscribable struct CollectLatestObservable{T, S, R, F} <: Subscribable{R}
    sources   :: S
    mappingFn :: F
end

function on_subscribe!(observable::CollectLatestObservable{L}, actor::A) where { L, A }
    sources = observable.sources
    storage = similar(sources, L)
    wrapper = CollectLatestObservableWrapper(L, actor, storage, observable.mappingFn)
    W       = typeof(wrapper)

    for index in CartesianIndices(axes(sources))
        @inbounds wrapper.subscriptions[index] = subscribe!(sources[index], CollectLatestObservableInnerActor{L, typeof(index), W}(index, wrapper))
        if cstatus(wrapper, index) === true && vstatus(wrapper, index) === false
            dispose(wrapper)
            break
        end
    end

    if all(wrapper.cstatus)
        dispose(wrapper)
    end

    return CollectLatestSubscription(wrapper)
end

##

struct CollectLatestSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: CollectLatestSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::CollectLatestSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::CollectLatestObservable{D}) where D  = print(io, "CollectLatestObservable($D)")
Base.show(io::IO, ::CollectLatestSubscription)           = print(io, "CollectLatestSubscription()")


text:
using PyPlot, SparseArrays, LinearAlgebra # Packages used

rows = [1,1,1,2,2,2,3,3,4,4,4,4,5,6,6,7,9,10,10,11,11,11]
cols = [2,10,6,3,1,11,1,2,5,3,2,7,6,7,5,6,10,1,9,1,4,9]
AT = sparse(cols, rows, 1, 11, 11)   # Use transpose, better for finding neighbors
spy(AT, marker=".", markersize=20);

function sparse_dfs(AT, start)
    visited = falses(size(AT,1))
    function visit(ivertex)
        visited[ivertex] = true
        println("Visiting vertex #$ivertex")
        for nb in findall(AT[:,ivertex] .!= 0)
            if !visited[nb]
                visit(nb)
            end
        end
    end
    visit(start)
    return nothing
end

sparse_dfs(AT, 1)


text:
include("parameters/parameters.jl")
"""
    AbstractItem
"""
abstract type AbstractItem end

"""
    Item <: AbstractItem

# Description

An immutable containing information about an item (a question in a test), e.g. `id::String` the item identifier, and the field `parameters::AbstractParameter` which accepts a mutable item parameter object

# Fields

  - **`idx::Int64`**: An integer that identifies the item in this session.
  - **`id::String`**: A string that identifies the examinee.
  - **`content::Vector{String}`**: A string vector containing the content features of an item.
  - **`parameters::AbstractParameters`**: A generic item parameters object.

# Factories
    Item(idx::Int64, id::String, content::Vector{String}, parameters::AbstractParameters) = new(idx, id, content, parameters)
    Item(idx::Int64, id::String, parameters::AbstractParameters) = new(idx, id, String[], parameters)
    Item(idx::Int64, parameters::AbstractParameters) = new(idx, string(idx), String[], parameters)

Creates a new generic item with custom index, id, (content features) and item parameters.

# Random initilizers
    Item(idx::Int64, id::String) = new(idx, id, String[], Parameters1PL())
    Item(idx::Int64, id::String, content::String) = new(idx, id, content, Parameters1PL())

Randomly generates a new generic calibrated item with custom index, id, content features and default 1PL item parameters
(Look at (`Parameters1PL`)[#Psychometrics.Parameters1PL] for the defaults).
"""
struct Item <: AbstractItem
    idx::Int64
    id::String
    content::Vector{String}
    parameters::AbstractParameters

    # Factories
    Item(idx::Int64, id::String, content::Vector{String}, parameters::AbstractParameters) = new(idx, id, content, parameters)
    Item(idx::Int64, id::String, parameters::AbstractParameters) = new(idx, id, String[], parameters)
    Item(idx::Int64, parameters::AbstractParameters) = new(idx, string(idx), String[], parameters)

    # Random default initilizers
    Item(idx::Int64, id::String) = new(idx, id, String[], Parameters1PL())
    Item(idx::Int64, id::String, content::String) = new(idx, id, content, Parameters1PL())
end

function Base.copy(item::Item)
    i = Item(item.idx, item.id, item.parameters)
    return i::Item
end

"""
    get_parameters(item::AbstractItem)
"""
function get_parameters(item::AbstractItem)
    item.parameters
end

"""
    get_item_by_id(item_id::String, items::Vector{<:AbstractItem})

# Description

It returns the item with index `item_id` from a vector of <:AbstractItem.
"""
function get_item_by_id(item_id::String, items::Vector{<:AbstractItem})
    filter(i -> i.id == item_id, items)[1]
end

"""
    get_parameters_vals(item::AbstractItem)
"""
function get_parameters_vals(item::AbstractItem)
    _get_parameters_vals(item.parameters)
end

"""
    empty_chain!(item::AbstractItem)
"""
function empty_chain!(item::AbstractItem)
    _empty_chain!(item.parameters)
end

"""
    chain_append!(item::AbstractItem)
"""
function chain_append!(item::AbstractItem)
    _chain_append!(item.parameters)
end


text:
"""
## rod_km!

This subroutine forms the stiffness matrix of a 1-d "rod" fin_el.

### Function
```julia
rod_km!(km, ea, length)
```

### Arguments
```julia
* km::Matrix{Float64}       : Element stiffness matrix (Updated)
* ea::Float64               : Element stiffness
* ell::Float64              : Element length
```
"""
function rod_km!(km::Matrix{Float64}, ea::Float64, ell::Float64)
  #
  # This subroutine forms the stiffness matrix of a 1-d "rod" fin_el.
  #
 km = [1.0 -1.0;-1.0 1.0]
 km *= ea/ell
end


text:
function fastnorm(u)
  # dest[1] = ...
  n = length(u)
  T = eltype(u)
  s = zero(T)
  @fastmath @inbounds @simd for i in 1:n
      s += u[i]^2
  end
  @fastmath @inbounds return sqrt(s)
end

function numericRoot(residFnc::Function, measurement, parameters, x0::Vector{Float64})
  # function is being deprecated
  return (nlsolve(   (res, X) -> residFnc(res, measurement, parameters, X), x0, inplace=true )).zero
end

function freshSamples!(ccwl::CommonConvWrapper, N::Int=1)
  ccwl.measurement = getSample(ccwl.usrfnc!, N)
  nothing
end

function shuffleXAltD(X::Vector{Float64}, Alt::Vector{Float64}, d::Int, p::Vector{Int})
  # n = length(X)
  Y = deepcopy(Alt)
  for i in 1:d
    Y[p[i]] = X[i]
  end
  return Y
end

"""
    $(SIGNATURES)

Shuffle incoming X into random positions in fr.Y.
Shuffled fr.Y will be placed back into fr.X[:,fr.gwp.particleidx] upon fr.gwp.usrfnc(x, res).
"""
function shuffleXAltD!(ccwl::CommonConvWrapper, X::Vector{Float64})
  # populate defaults from existing values
  for i in 1:ccwl.xDim
    ccwl.cpt[Threads.threadid()].Y[i] = ccwl.cpt[Threads.threadid()].X[i, ccwl.cpt[Threads.threadid()].particleidx]
  end
  # populate as many measurment dimensions randomly for calculation
  for i in 1:ccwl.zDim
    ccwl.cpt[Threads.threadid()].Y[ccwl.cpt[Threads.threadid()].p[i]] = X[i]
  end
  nothing
end

function (ccw::CommonConvWrapper)(res::Vector{Float64}, x::Vector{Float64})
  shuffleXAltD!(ccw, x)
  ccw.params[ccw.varidx][:, ccw.cpt[Threads.threadid()].particleidx] = ccw.cpt[Threads.threadid()].Y
  ret = ccw.usrfnc!(res, ccw.cpt[Threads.threadid()].factormetadata, ccw.cpt[Threads.threadid()].particleidx, ccw.measurement, ccw.params[ccw.cpt[Threads.threadid()].activehypo]...) # optmize the view here, re-use the same memory
  return ret
end

function (ccw::CommonConvWrapper)(x::Vector{Float64})
  ccw.params[ccw.varidx][:, ccw.cpt[Threads.threadid()].particleidx] = x #ccw.Y
  ccw.usrfnc!(ccw.cpt[Threads.threadid()].res, ccw.cpt[Threads.threadid()].factormetadata, ccw.cpt[Threads.threadid()].particleidx, ccw.measurement, ccw.params[ccw.cpt[Threads.threadid()].activehypo]...)
end

function numericRootGenericRandomizedFnc!(
            ccwl::CommonConvWrapper{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwiseMinimize}
  #
  fill!(ccwl.cpt[Threads.threadid()].res, 0.0) # 1:frl.xDim
  r = optimize( ccwl, ccwl.cpt[Threads.threadid()].X[:, ccwl.cpt[Threads.threadid()].particleidx] ) # ccw.gg
  # TODO -- clearly lots of optmization to be done here
  ccwl.cpt[Threads.threadid()].Y[:] = r.minimizer
  ccwl.cpt[Threads.threadid()].X[:,ccwl.cpt[Threads.threadid()].particleidx] = ccwl.cpt[Threads.threadid()].Y
  nothing
end

"""
    $(SIGNATURES)

Solve free variable x by root finding residual function fgr.usrfnc(x, res)
randomly shuffle x dimensions if underconstrained by measurement z dimensions
small random perturbation used to prevent trivial solver cases, div by 0 etc.
result stored in fgr.Y
ccw.X must be set to memory ref the param[varidx] being solved, at creation of ccw
"""
function numericRootGenericRandomizedFnc!(
            ccwl::CommonConvWrapper{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwise}
  #
  ## TODO desperately needs cleaning up and refactoring
  # ststr = "thrid=$(Threads.threadid()), zDim=$(ccwl.zDim), xDim=$(ccwl.xDim)\n"
  # ccall(:jl_, Nothing, (Any,), ststr)
  if ccwl.zDim < ccwl.xDim && !ccwl.partial || testshuffle
    # less measurement dimensions than variable dimensions -- i.e. shuffle
    shuffle!(ccwl.cpt[Threads.threadid()].p)
    for i in 1:ccwl.xDim
      ccwl.cpt[Threads.threadid()].perturb[1:ccwl.zDim] = perturb*randn(ccwl.zDim)
      ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx] += ccwl.cpt[Threads.threadid()].perturb
      r = nlsolve(  ccwl,
                    ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx], # this is x0
                    inplace=true
                 )
      if r.f_converged
        shuffleXAltD!( ccwl, r.zero )
        break;
      else
        # TODO -- report on this bottleneck, useful for optimization of code
        # @show i, ccwl.p, ccwl.xDim, ccwl.zDim
        temp = ccwl.cpt[Threads.threadid()].p[end]
        ccwl.cpt[Threads.threadid()].p[2:end] = ccwl.cpt[Threads.threadid()].p[1:(end-1)]
        ccwl.cpt[Threads.threadid()].p[1] = temp
        if i == ccwl.xDim
          error("numericRootGenericRandomizedFnc could not converge, i=$(i), ccwl.usrfnc!=$(typeof(ccwl.usrfnc!))")
        end
      end
    end
    #shuffleXAltD!( ccwl, r.zero ) # moved up
  elseif ccwl.zDim >= ccwl.xDim && !ccwl.partial
    # equal or more measurement dimensions than variable dimensions -- i.e. don't shuffle
    ccwl.cpt[Threads.threadid()].perturb[1:ccwl.xDim] = perturb*randn(ccwl.xDim)
    ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim, ccwl.cpt[Threads.threadid()].particleidx] += ccwl.cpt[Threads.threadid()].perturb[1:ccwl.xDim] # moved up

    # str = "nlsolve, thrid_=$(Threads.threadid()), partidx=$(ccwl.cpt[Threads.threadid()].particleidx), X=$(ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim,ccwl.cpt[Threads.threadid()].particleidx])"
    # ccall(:jl_, Nothing, (Any,), str)
    r = nlsolve( ccwl, ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim,ccwl.cpt[Threads.threadid()].particleidx], inplace=true )
    # ccall(:jl_, Nothing, (Any,), "nlsolve.zero=$(r.zero)")
    if sum(isnan.(( r ).zero)) == 0
      ccwl.cpt[Threads.threadid()].Y[1:ccwl.xDim] = ( r ).zero
    else
      # TODO print this output as needed
      str = "ccw.thrid_=$(Threads.threadid()), got NaN, ccwl.cpt[Threads.threadid()].particleidx = $(ccwl.cpt[Threads.threadid()].particleidx), r=$(r)\n"
      @info str
      ccall(:jl_, Nothing, (Any,), str)
      ccall(:jl_, Nothing, (Any,), ccwl.usrfnc!)
      for thatlen in 1:length(ccwl.params)
        str = "thatlen=$thatlen, ccwl.params[thatlen][:, ccwl.cpt[Threads.threadid()].particleidx]=$(ccwl.params[thatlen][:, ccwl.cpt[Threads.threadid()].particleidx])\n"
        ccall(:jl_, Nothing, (Any,), str)
        @warn str
      end
    end
  elseif ccwl.partial
    # improve memory management in this function
    ccwl.cpt[Threads.threadid()].p[1:length(ccwl.usrfnc!.partial)] = Int[ccwl.usrfnc!.partial...] # TODO -- move this line up and out of inner loop
    r = nlsolve(  ccwl,
                  ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx], # this is x0
                  inplace=true
               )
    shuffleXAltD!( ccwl, r.zero )
  else
    ccall(:jl_, Nothing, (Any,), "ERROR: Unresolved numeric solve case")
    error("Unresolved numeric solve case")
  end
  ccwl.cpt[Threads.threadid()].X[:,ccwl.cpt[Threads.threadid()].particleidx] = ccwl.cpt[Threads.threadid()].Y
  nothing
end

"""
    $(SIGNATURES)

Set variable(s) `sym` of factor graph to be marginalized -- i.e. not be updated by inference computation.
"""
function setfreeze!(dfg::AbstractDFG, sym::Symbol)
  if !isInitialized(dfg, sym)
    @warn "Vertex $(sym) is not initialized, and won't be frozen at this time."
    return nothing
  end
  vert = DFG.getVariable(dfg, sym)
  data = solverData(vert)
  data.ismargin = true
  nothing
end
function setfreeze!(dfg::AbstractDFG, syms::Vector{Symbol})
  for sym in syms
    setfreeze!(dfg, sym)
  end
end

"""
    $(SIGNATURES)

Freeze nodes that are older than the quasi fixed-lag length defined by `fg.qfl`, according to `fg.fifo` ordering.

Future:
- Allow different freezing strategies beyond fifo.
"""
function fifoFreeze!(dfg::G)::Nothing where G <: AbstractDFG
  if DFG.getSolverParams(dfg).qfl == 0
    @warn "Quasi fixed-lag is enabled but QFL horizon is zero. Please set a valid window with FactoGraph.qfl"
  end

  # the fifo history
  tofreeze = DFG.getAddHistory(dfg)[1:(end-DFG.getSolverParams(dfg).qfl)]
  if length(tofreeze) == 0
      @info "[fifoFreeze] QFL - no nodes to freeze."
      return nothing
  end
  @info "[fifoFreeze] QFL - Freezing nodes $(tofreeze[1]) -> $(tofreeze[end])."
  setfreeze!(dfg, tofreeze)
  nothing
end

"""
    $(SIGNATURES)

Return all factors currently registered in the workspace.
"""
function getCurrentWorkspaceFactors()::Vector{Type}
    return [
        subtypes(IncrementalInference.FunctorSingleton)...,
        subtypes(IncrementalInference.FunctorPairwise)...,
        subtypes(IncrementalInference.FunctorPairwiseMinimize)...];
end

"""
    $(SIGNATURES)

Return all variables currently registered in the workspace.
"""
function getCurrentWorkspaceVariables()::Vector{Type}
    return subtypes(IncrementalInference.InferenceVariable);
end

"""
    $SIGNATURES

Reset the state of all variables in a clique to not initialized.

Notes
- resets numberical values to zeros.

Dev Notes
- TODO not all kde manifolds will initialize to zero.
"""
function resetCliqSolve!(dfg::G,
                         treel::BayesTree,
                         cliq::Graphs.ExVertex;
                         solveKey::Symbol=:default)::Nothing where G <: AbstractDFG
  #
  cda = getData(cliq)
  vars = getCliqVarIdsAll(cliq)
  for varis in vars
    resetVariable!(dfg, varis, solveKey=solveKey)
  end
  prnt = getParent(treel, cliq)
  if length(prnt) > 0
    setCliqUpInitMsgs!(prnt[1], cliq.index, TempBeliefMsg())
  end
  cda.upMsg = Dict{Symbol, BallTreeDensity}()
  cda.dwnMsg = Dict{Symbol, BallTreeDensity}()
  cda.upInitMsgs = Dict{Int, TempBeliefMsg}()
  cda.downInitMsg = TempBeliefMsg()
  setCliqStatus!(cliq, :null)
  setCliqDrawColor(cliq, "")
  return nothing
end

function resetCliqSolve!(dfg::G,
                         treel::BayesTree,
                         frt::Symbol;
                         solveKey::Symbol=:default  )::Nothing where G <: AbstractDFG
  #
  resetCliqSolve!(dfg, treel, getCliq(treel, frt), solveKey=solveKey)
end

"""
    $SIGNATURES

Inverse solve of predicted noise value and returns the associated "measured" noise value (also used as starting point for the solve).
"""
function solveFactorMeasurements(dfg::AbstractDFG,
                                 fctsym::Symbol  )
  #
  fcto = getFactor(dfg, fctsym)
  varsyms = fcto._variableOrderSymbols
  vars = map(x->getPoints(getKDE(dfg,x)), varsyms)
  fcttype = getFactorType(fcto)
  zDim = getData(fcto).fnc.zDim

  N = size(vars[1])[2]
  res = zeros(zDim)
  ud = FactorMetadata()
  meas = getSample(fcttype, N)
  meas0 = deepcopy(meas[1])

  function makemeas!(i, meas, dm)
    meas[1][:,i] = dm
    return meas
  end

  ggo = (i, dm) -> fcttype(res,ud,i,makemeas!(i, meas, dm),vars...)
  # ggo(1, [0.0;0.0])

  for idx in 1:N
    retry = 10
    while 0 < retry
      if isa(fcttype, FunctorPairwiseMinimize)
        r = optimize((x) -> ggo(idx,x), meas[1][:,idx]) # zeros(zDim)
        retry -= 1
        if !r.g_converged
          nsm = getSample(fcttype, 1)
          for count in 1:length(meas)
            meas[count][:,idx] = nsm[count][:,idx]
          end
        else
          break
        end
      elseif isa(fcttype, FunctorPairwise)
        ggnl = (rs, dm) -> fcttype(rs,ud,idx,makemeas!(idx, meas, dm),vars...)
        r = nlsolve(ggnl, meas[1][:,idx])
        break
      elseif isa(fcttype, FunctorSingleton)
        # assuming no partials at this point
        meas[1][:,:] .= vars[1][:,:]
        break
      end
    end
    # @assert meas[1][:,idx] == r.minimizer
  end

  # Gadfly.plot(z=(x,y)->ggo(1,[x;y]), xmin=[-pi],xmax=[pi],ymin=[-100.0],ymax=[100.0], Geom.contour)
  return meas[1], meas0
end

#


text:
include("reproducibility.jl")

using SymPy
using QSWalk
using SparseArrays
using LinearAlgebra

x = symbols("x")
w = symbols("w")

for n = 3:3
    # 3 - nodes LQSW example
    a = spzeros(n, n)
    a[1,2:n] .= 1
    @show n
    l = local_lind(a)
    h = a + transpose(a) # a is upper-triangular so this is enough
    s_h = evolve_generator(h, [zeros(size(h)...)])
    s_l = evolve_generator(zeros(size(h)...), l)
    s_h = Matrix(s_h)
    s_l = Matrix(s_l)
    #println(s_h)

    #println(s_l)

    #continue
    s = (1-w)*s_h+w*s_l
    det_matrix = s - Diagonal(ones(n^2))*x
    det_s = det_matrix.det()
    det_s_der = diff(det_s, x)
    polynomial = det_s_der.subs(x, 0)
    @show polynomial
    @show solve(polynomial, w)
end


text:
module ClimaCoreMakie

import Makie: Makie, @recipe
import ClimaCore

@recipe(Viz, object) do scene
    Makie.Attributes(; colormap = :balance)
end

function level_mesh(field::ClimaCore.Fields.SpectralElementField2D, level = 1)
    space = axes(field)
    Ni, Nj, _, _, Nh = size(space.local_geometry)
    @assert Ni == Nj

    cart_coords =
        ClimaCore.Geometry.Cartesian123Point.(
            ClimaCore.Fields.coordinate_field(space),
            Ref(space.global_geometry),
        )

    triangle_count = Nh * (Nj - 1) * (Ni - 1) * 2
    FaceType = Makie.GeometryBasics.GLTriangleFace
    PointType = Makie.GeometryBasics.Point3f

    faces = Array{FaceType}(undef, triangle_count)
    vertices = Array{PointType}(undef, triangle_count * 3)

    x1_data = ClimaCore.Fields.todata(cart_coords.x1)
    x2_data = ClimaCore.Fields.todata(cart_coords.x2)
    x3_data = ClimaCore.Fields.todata(cart_coords.x3)

    L = LinearIndices((1:Ni, 1:Nj))
    idx = 0
    for h in 1:Nh
        x1_slab = ClimaCore.slab(x1_data, level, h)
        x2_slab = ClimaCore.slab(x2_data, level, h)
        x3_slab = ClimaCore.slab(x3_data, level, h)
        for j in 1:(Nj - 1), i in 1:(Ni - 1), t in 1:2
            if t == 1
                # quad triangle 1
                I = L[i, j] + Ni * Nj * (h - 1)
                J = L[i + 1, j] + Ni * Nj * (h - 1)
                K = L[i, j + 1] + Ni * Nj * (h - 1)
                faces[idx + 1] = Makie.GeometryBasics.GLTriangleFace(I, J, K)
                vertices[idx * 3 + 1] = Makie.GeometryBasics.Point3f(
                    x1_slab[i, j],
                    x2_slab[i, j],
                    x3_slab[i, j],
                )
                vertices[idx * 3 + 2] = Makie.GeometryBasics.Point3f(
                    x1_slab[i + 1, j],
                    x2_slab[i + 1, j],
                    x3_slab[i + 1, j],
                )
                vertices[idx * 3 + 3] = Makie.GeometryBasics.Point3f(
                    x1_slab[i, j + 1],
                    x2_slab[i, j + 1],
                    x3_slab[i, j + 1],
                )
            else
                # quad triangle 2
                I = L[i + 1, j] + Ni * Nj * (h - 1)
                J = L[i + 1, j + 1] + Ni * Nj * (h - 1)
                K = L[i, j + 1] + Ni * Nj * (h - 1)
                faces[idx + 1] = Makie.GeometryBasics.GLTriangleFace(I, J, K)
                vertices[idx * 3 + 1] = Makie.GeometryBasics.Point3f(
                    x1_slab[i + 1, j],
                    x2_slab[i + 1, j],
                    x3_slab[i + 1, j],
                )
                vertices[idx * 3 + 2] = Makie.GeometryBasics.Point3f(
                    x1_slab[i + 1, j + 1],
                    x2_slab[i + 1, j + 1],
                    x3_slab[i + 1, j + 1],
                )
                vertices[idx * 3 + 3] = Makie.GeometryBasics.Point3f(
                    x1_slab[i, j + 1],
                    x2_slab[i, j + 1],
                    x3_slab[i, j + 1],
                )
            end
            idx += 1
        end
    end
    return Makie.GeometryBasics.Mesh(vertices, faces)
end

Makie.plottype(::ClimaCore.Fields.SpectralElementField2D) =
    Viz{<:Tuple{ClimaCore.Fields.SpectralElementField2D}}

function Makie.plot!(
    plot::Viz{<:Tuple{ClimaCore.Fields.SpectralElementField2D}},
)
    # retrieve the field to plot
    field = plot[:object][]

    if !(eltype(field) <: Union{Float64, Float32})
        error("plotting only implemented for F64, F32 scalar fields")
    end
    space = ClimaCore.axes(field)

    cart_coords =
        ClimaCore.Geometry.Cartesian123Point.(
            ClimaCore.Fields.coordinate_field(space),
            Ref(space.global_geometry),
        )
    vertices = hcat(
        vec(parent(cart_coords.x1)),
        vec(parent(cart_coords.x2)),
        vec(parent(cart_coords.x3)),
    )
    triangles = hcat(ClimaCore.Spaces.triangulate(space)...)

    colors = vec(parent(field))
    Makie.mesh!(
        plot,
        vertices,
        triangles,
        color = colors,
        colormap = plot[:colormap],
        shading = false,
    )
end

end # module


text:
"""
    r2(model::LinearMixedModel; conditional=false)
    r²(model::LinearMixedModel; conditional=false)

Coefficient of determination (R-squared).

R² is very non trivial for mixed models for a number of reasons
and the measures here are particularly simplistic. For `conditional=true`,
the Pearson correlation between the fitted and observed values is simply
squared, in line with the relationship between the coefficient of determination
and the Pearson correlation for classical OLS models. For `conditional=false`,
new predicted values are generated based on only the fixed-effects estimates
and the correlation between these predictions and the observed values is again
squared.

In both cases, it is important to note that despite the usual naming convention
(R² for the coefficient of determination and r for Pearson's correlation
coefficient), these quantitaties are not defined in terms of the other.
Instead, the correlation is a standardized measure of covariance and the
coefficient of determination is measured relative to a null model. (The total
sum of squares is in some sense the squared residual error of the null model.)
Even generalizations of the coefficient of determination that define this value
in terms of the likelihoods of an intercept-only model (i.e. the sample mean)
and the likelihood of the fitted model fail to generalize to the mixed-model
case: should "intercept-only" mean truly just the fixed effect (i.e. comparing
the likelihoods of a classical OLS and a mixed-effects model)? or should it
also include an random intercept for each grouping variable? What is the null
model in this case, when determining whether each grouping variable contributes
to overall fit? There is no single, clear, agreed-upon answer for these
questions.  The bottom line is that there are many potential ways to define a
coefficient of determination for linear mixed models (and even more for the
generalzed case which combines all the problems of R² for GLM and for LMM) and
none of them have all the properties of the coefficient of determination for
classical OLS regression.

There are more advanced approximations (see e.g. Nakagawa and colleagues' 2013
2017 papers) but the fundamental philosophical issues above remain.
In the future, additional approximations may be supported, but this is not a
high priority.  Pull requests are welcome.

For more information, see the [GLMM FAQ](https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#how-do-i-compute-a-coefficient-of-determination-r2-or-an-analogue-for-glmms)
"""
function StatsBase.r2(obj::LinearMixedModel; conditional=false)
    # The conditional case is approximately equivalent to
    # 1 - var(residuals(obj)) / var(obj.y)
    ŷ = conditional ? fitted(obj) : obj.X * obj.β
    return cor(ŷ, obj.y)^2
end

"""
    adjr2(model::LinearMixedModel; conditional=false)
    adjr²(model::LinearMixedModel; conditional=false)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as ``(1 - (1-R^2)(n-1)/(n-p))``, with ``R^2``
the coefficient of determination, ``n`` the number of observations, and ``p`` the number of
coefficients (including the intercept). This definition is generally known as the Wherry Formula I.

For mixed-models, we use this same adjustment for the R² value computed with
[`r2`](@ref). `conditional` whether or not the random effects are taken into
account when computing the quality of the model fit. They are nonetheless
included in the adjustment (i.e. as number ``p`` of parameters).

R² is not a clearly defined concept for mixed models, nor does it have all the
properties typically expected of the coefficient of determination. See
[`r2`](@ref) for more information.
"""
function StatsBase.adjr2(model::LinearMixedModel; conditional=false)
    R² = r2(model; conditional=conditional)
    n = nobs(model)
    p = dof(model)
    return 1 - (1 - R²) * (n - 1) / (n - p)
end


text:
println("Launching")
using MocosSimLauncher

launch(ARGS)


text:
#=
  From BLOG example/bayesian-linear-regression.blog
  """
  A model for Bayesian linear regression

  @author yiwu
  @since 2015-07-02
  """

  Note: This is 0-based i.e.
     x y
     0 0
     1 1
     2 4
     3 6
     4 ?

  Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

           m    1.9616    0.4171     0.0042    0.0218   198.2693    1.0229       25.6327
           b   -0.1645    0.7370     0.0074    0.0321   351.3579    1.0083       45.4244
       sigma    0.9349    0.4702     0.0047    0.0309    86.5685    1.0103       11.1918
          y4    7.6895    1.6275     0.0163    0.0742   229.2941    1.0250       29.6437
     post[1]   -0.1513    1.2702     0.0127    0.0441   765.0960    1.0014       98.9135
     post[2]    1.8138    1.1771     0.0118    0.0365   910.5267    1.0006      117.7152
     post[3]    3.7861    1.1241     0.0112    0.0386   667.4734    1.0052       86.2926
     post[4]    5.7280    1.3455     0.0135    0.0475   517.3522    1.0120       66.8846
     post[5]    7.7025    1.6070     0.0161    0.0705   311.8366    1.0193       40.3150
     post[6]    9.6650    1.9211     0.0192    0.0894   276.0979    1.0201       35.6946
     post[7]   11.5449    2.1770     0.0218    0.1019   276.8155    1.0169       35.7874
     post[8]   13.5644    2.6050     0.0260    0.1283   236.7469    1.0168       30.6072
     post[9]   15.5463    2.9901     0.0299    0.1471   210.4787    1.0206       27.2112
    post[10]   17.4897    3.3944     0.0339    0.1688   230.5202    1.0202       29.8022
    post[11]   19.4406    3.7283     0.0373    0.1901   214.3265    1.0238       27.7087

   Cf ~/webppl/bayesian_linear_regression.wppl

=#

using Turing
include("jl_utils.jl")

@model function bayesian_linear_regression(y=[0,1,4,6])
    m ~ Normal(0,2)
    b ~ Normal(0,2)
    sigma ~ Gamma(1,1)
    for i in 0:length(y)-1
        y[i+1] ~ Normal(m * i + b, sigma)
    end

    # Posterior
    y4 ~ Normal(m * (4) + b, sigma)
    post = Vector{Real}(undef, n)
    for i in 0:10
        post[i+1] ~ Normal(m * i + b, sigma)
    end
end

y = [0,1,4,6]
model = bayesian_linear_regression(y)

num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
# chs = sample(model, IS(), 10_000)
#chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)

# chs = sample(model, SGLD(), 10_000)
chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank