text:
module Shapes

using IterTools: partition
using LinearAlgebra
using StaticArrays: SVector, SMatrix

export AbstractShape, AbstractLineShape,
    HalfPlane, HalfSpace, Point, Polygon, Polyline, Vertex,
    circle, edges, polygon, rotate, translate, vertices

const flatten = Iterators.flatten

# """
#     Vector{D,T}

# A geometric vector in `D` dimensions of data type `T`.

# """
# struct Vector{D,T}
#     data::NTuple{D,T}
# end
# Vector(coordinates...) = Vector(promote(coordinates...))
# Vector(n::T) where {T<:Number} = Vector{1,T}((n,))
# Vector(iterable) = Vector(iterable...)
# Base.convert(::Type{Vector{D,T}}, p::Vector{D,S}) where {D,T,S} = Vector(convert(NTuple{D,T}, p.data))
# Base.promote_rule(::Type{Vector{D,T}}, ::Type{Vector{D,S}}) where {D,T,S} = Vector{D,promote_type(T, S)}

# @inline Base.getproperty(a::Vector, v::Symbol) = _getproperty(a, Val(v))
# @inline _getproperty(a::Vector, ::Val{v}) where {v} = Core.getfield(a, v)
# @inline _getproperty(a::Vector, ::Val{:x}) = Core.getfield(a, :data)[1]
# @inline _getproperty(a::Vector, ::Val{:y}) = Core.getfield(a, :data)[2]
# @inline _getproperty(a::Vector, ::Val{:z}) = Core.getfield(a, :data)[3]
# Base.iterate(a::Vector) = iterate(Core.getfield(a, :data))
# Base.iterate(a::Vector, n) = iterate(Core.getfield(a, :data), n)
# @inline Base.getindex(a::Vector, idx) = Base.getindex(Core.getfield(a, :data), idx)
# Base.length(a::Vector) = length(Core.getfield(a, :data))
# Base.lastindex(a::Vector) = lastindex(Core.getfield(a, :data))

##################################################
# Basic shapes
##################################################

abstract type AbstractShape{D,T} end

"""
    AbstractPrimitiveShape{D,T}

A primitive shape in `D` dimensions of data type `T`.

"""
abstract type AbstractPrimitiveShape{D,T} <: AbstractShape{D,T} end

"""
    AbstractLineShape{N,D,T}

A shape consisting of `N` vertices connected by line segments in `D`
dimensions of data type `T`.

"""
abstract type AbstractLineShape{N,D,T} <: AbstractShape{D,T} end

"""
    Vertex{D,T}

A vertex in `D` dimensions of data type `T`.

"""
struct Vertex{D,T} <: AbstractPrimitiveShape{D,T}
    data::NTuple{D,T}
end
Vertex(coordinates...) = Vertex(promote(coordinates...))
Vertex(n::T) where {T<:Number} = Vertex{1,T}((n,))
Vertex(iterable) = Vertex(iterable...)
Base.convert(::Type{Vertex{D,T}}, p::Vertex{D,S}) where {D,T,S} = Vertex(convert(NTuple{D,T}, p.data))
Base.promote_rule(::Type{Vertex{D,T}}, ::Type{Vertex{D,S}}) where {D,T,S} = Vertex{D,promote_type(T, S)}

@inline Base.getproperty(a::Vertex, v::Symbol) = _getproperty(a, Val(v))
@inline _getproperty(a::Vertex, ::Val{v}) where {v} = Core.getfield(a, v)
@inline _getproperty(a::Vertex, ::Val{:x}) = Core.getfield(a, :data)[1]
@inline _getproperty(a::Vertex, ::Val{:y}) = Core.getfield(a, :data)[2]
@inline _getproperty(a::Vertex, ::Val{:z}) = Core.getfield(a, :data)[3]
Base.iterate(a::Vertex) = iterate(Core.getfield(a, :data))
Base.iterate(a::Vertex, n) = iterate(Core.getfield(a, :data), n)
@inline Base.getindex(a::Vertex, idx) = Base.getindex(Core.getfield(a, :data), idx)
Base.length(a::Vertex) = length(Core.getfield(a, :data))
Base.lastindex(a::Vertex) = lastindex(Core.getfield(a, :data))

"""
    Polyline{N,D,T}

A collection of line segments connecting `N` vertices, in `D`
dimensions of data type `T`.

"""
struct Polyline{N,D,T} <: AbstractLineShape{N,D,T}
    data::NTuple{N,Vertex{D,T}}
end
Polyline(iterable) = Polyline(Tuple(iterable))
Polyline(vertices::Vertex{D,T}...) where {D,T} = Polyline(vertices)

const Point{D,T} = Polyline{1,D,T}
const LineSegment{D,T} = Polyline{2,D,T}

"""
    Polygon{N,D,T}

Similar to `Polyline`, except the last and first vertices are also
connected by a line segment.

"""
struct Polygon{N,D,T} <: AbstractLineShape{N,D,T}
    data::NTuple{N,Vertex{D,T}}
end
Polygon(iterable) = Polygon(Tuple(iterable))
Polygon(vertices::Vertex{D,T}...) where {D,T} = Polygon(vertices)

Base.iterate(a::AbstractLineShape) = iterate(a.data)
Base.iterate(a::AbstractLineShape, n) = iterate(a.data, n)
Base.getindex(a::AbstractLineShape, idx) = getindex(a.data, idx)
Base.length(a::AbstractLineShape) = length(a.data)
Base.lastindex(a::AbstractLineShape) = lastindex(a.data)

circle(c::Vertex{2}, r::Real; samples=20, θ=0.) = Polygon(Vertex(c.x + r*cosd(360*i/samples + θ), c.y + r*sind(360*i/samples + θ)) for i in 0:samples-1)

"""
    HalfSpace{D,T}

A shape that separates a space of `D` dimensions into two.

"""
struct HalfSpace{D,T} <: AbstractShape{D,T}
    vertex::Vertex{D,T}
    normal::SVector{D,T}
end

const HalfLine{T} = HalfSpace{1,T}
const HalfPlane{T} = HalfSpace{2,T}

##################################################
# BoundingBox
##################################################

function axes(::AbstractShape{D,T}) where {D,T}
    m = SMatrix{D,D}(Matrix(one(T)*I, D, D))
    return [m[:,i] for i in 1:D]
end

##################################################
# BoundingBox
##################################################

"""
    BoundingBox{D,T}

A bounding box in `D` dimensions of data type `T`.

"""
struct BoundingBox{D,T}
    data::NTuple{D,Pair{T,T}}
end
# BoundingBox(iterable) = BoundingBox(Tuple(iterable))
BoundingBox(bounds::Pair{T,T}...) where T = BoundingBox(bounds)

@inline Base.getproperty(a::BoundingBox, v::Symbol) = _getproperty(a, Val(v))
@inline _getproperty(a::BoundingBox, ::Val{v}) where v = Core.getfield(a, v)
@inline _getproperty(a::BoundingBox, ::Val{:x}) = Core.getfield(a, :data)[1]
@inline _getproperty(a::BoundingBox, ::Val{:y}) = Core.getfield(a, :data)[2]
@inline _getproperty(a::BoundingBox, ::Val{:z}) = Core.getfield(a, :data)[3]
Base.iterate(a::BoundingBox) = iterate(Core.getfield(a, :data))
Base.iterate(a::BoundingBox, n) = iterate(Core.getfield(a, :data), n)
Base.getindex(a::BoundingBox, idx) = getindex(Core.getfield(a, :data), idx)
Base.length(a::BoundingBox) = length(Core.getfield(a, :data))
Base.lastindex(a::BoundingBox) = lastindex(a.data)

BoundingBox(b::BoundingBox) = b
BoundingBox(a::BoundingBox{1}, b::BoundingBox{1}) = BoundingBox(
    min(first(a.x), first(b.x)) => max(last(a.x), last(b.x)))
BoundingBox(a::BoundingBox{2}, b::BoundingBox{2}) = BoundingBox(
    min(first(a.x), first(b.x)) => max(last(a.x), last(b.x)),
    min(first(a.y), first(b.y)) => max(last(a.y), last(b.y)))
BoundingBox(a::BoundingBox{3}, b::BoundingBox{3}) = BoundingBox(
    min(first(a.x), first(b.x)) => max(last(a.x), last(b.x)),
    min(first(a.y), first(b.y)) => max(last(a.y), last(b.y)),
    min(first(a.z), first(b.z)) => max(last(a.z), last(b.z)))
BoundingBox(a::BoundingBox{D}, b::BoundingBox{D}) where {D} = BoundingBox(
    min(first(ax), first(bx)) => max(last(ax), last(bx)) for (ax,bx) in zip(a, b))

BoundingBox(a::Vertex{2}) = BoundingBox(a.x => a.x, a.y => a.y)
BoundingBox(a::Vertex{3}) = BoundingBox(a.x => a.x, a.y => a.y, a.z => a.z)
BoundingBox(a::Vertex) = BoundingBox(x => x for x in a)

function BoundingBox(a::AbstractLineShape{N,1}) where {N}
    xmin, xmax = Inf, -Inf
    for v in vertices(a)
        v.x < xmin && (xmin = v.x)
        v.x > xmax && (xmax = v.x)
    end
    return BoundingBox(xmin => xmax)
end

function BoundingBox(a::AbstractLineShape{N,2}) where {N}
    xmin, xmax = Inf, -Inf
    ymin, ymax = Inf, -Inf
    for v in vertices(a)
        v.x < xmin && (xmin = v.x)
        v.x > xmax && (xmax = v.x)
        v.y < ymin && (ymin = v.y)
        v.y > ymax && (ymax = v.y)
    end
    return BoundingBox(xmin => xmax, ymin => ymax)
end

function BoundingBox(a::AbstractLineShape{N,3}) where {N}
    xmin, xmax = Inf, -Inf
    ymin, ymax = Inf, -Inf
    zmin, zmax = Inf, -Inf
    for v in vertices(a)
        v.x < xmin && (xmin = v.x)
        v.x > xmax && (xmax = v.x)
        v.y < ymin && (ymin = v.y)
        v.y > ymax && (ymax = v.y)
        v.z < zmin && (zmin = v.z)
        v.z > zmax && (zmax = v.z)
    end
    return BoundingBox(xmin => xmax, ymin => ymax, zmin => zmax)
end

function BoundingBox(a::AbstractLineShape{N,D}) where {N,D}
    mins = fill(Inf, D)
    maxs = fill(-Inf, D)
    for v in vertices(a), d in 1:D
        v[d] < mins[d] && (mins[d] = v[d])
        v[d] > maxs[d] && (maxs[d] = v[d])
    end
    return BoundingBox(map(Pair, mins, maxs)...)
end

function BoundingBox(s::HalfSpace{D,T}) where {D,T}
    mins = fill(Inf, D)
    maxs = fill(-Inf, D)
    axs = axes(s)
    for i in 1:D
        a = axs[i]
        if iszero(a × s.normal)
            x = a ⋅ vector(s.vertex)
            if a ⋅ s.normal >= 0
                mins[i] = -Inf
                maxs[i] = x
                continue
            end
            mins[i] = x
            maxs[i] = Inf
            continue
        end
        mins[i] = -Inf
        maxs[i] = Inf
    end
    return BoundingBox(map(Pair, mins, maxs)...)
end

intersects(a::BoundingBox{1}, b::BoundingBox{1}) =
    last(a.x) > first(b.x) && last(b.x) > first(a.x)
intersects(a::BoundingBox{2}, b::BoundingBox{2}) =
    last(a.x) > first(b.x) && last(b.x) > first(a.x) &&
    last(a.y) > first(b.y) && last(b.y) > first(a.y)
intersects(a::BoundingBox{3}, b::BoundingBox{3}) =
    last(a.x) > first(b.x) && last(b.x) > first(a.x) &&
    last(a.y) > first(b.y) && last(b.y) > first(a.y) &&
    last(a.z) > first(b.z) && last(b.z) > first(a.z)
intersects(a::BoundingBox{D}, b::BoundingBox{D}) where {D} = all(
    last(ax) > first(bx) && last(bx) > first(ax) for (ax,bx) in zip(a, b))

##################################################
# Basic algebra
##################################################

Base.isapprox(a::Vertex{1}, b::Vertex{1}) = isapprox(a.x, b.x)
Base.isapprox(a::Vertex{2}, b::Vertex{2}) = isapprox(a.x, b.x) && isapprox(a.y, b.y)
Base.isapprox(a::Vertex{3}, b::Vertex{3}) = isapprox(a.x, b.x) && isapprox(a.y, b.y) && isapprox(a.z, b.z)
Base.isapprox(a::Vertex{D}, b::Vertex{D}) where {D} = all(isapprox(ax, bx) for (ax,bx) in zip(a, b))

"""
    vector(a::Vertex)

Return a vector pointing from the origin to `a`.

"""
vector(a::Vertex{1}) = SVector(a.x)
vector(a::Vertex{2}) = SVector(a.x, a.y)
vector(a::Vertex{3}) = SVector(a.x, a.y, a.z)
vector(a::Vertex) = SVector(a.data...)

"""
    vector(a::Vertex, b::Vertex)

Return a vector pointing from `a` to `b`.

"""
vector(a::Vertex{1}, b::Vertex{1}) = SVector(b.x - a.x)
vector(a::Vertex{2}, b::Vertex{2}) = SVector(b.x - a.x, b.y - a.y)
vector(a::Vertex{3}, b::Vertex{3}) = SVector(b.x - a.x, b.y - a.y, b.z - a.z)
vector(a::Vertex{N}, b::Vertex{N}) where N = SVector(map(-, b, a)...)

"""
    vector(p::Point)

Return a vector pointing from the origin to `p`.

"""
vector(p::Point) = vector(first(p))

"""
    vector(l::LineSegment)

Return a vector pointing along the line segment `l`.

"""
vector(l::LineSegment) = vector(first(l), last(l))

"""
    vertices(shape)
"""
vertices(p::Vertex) = (p,)
vertices(l::AbstractLineShape) = l.data

"""
    edges(shape)
"""
edges(p::Vertex{D,T}) where {D,T} = Polyline{2,D,T}[]

function edges(l::Polyline{N,D,T}) where {N,D,T}
    result = Base.Vector{Polyline{2,D,T}}(undef, N - 1)
    for i in 1:N-1
        result[i] = Polyline(l[i], l[i+1])
    end
    return result
end

function edges(l::Polygon{N,D,T}) where {N,D,T}
    result = Base.Vector{Polyline{2,D,T}}(undef, N)
    for i in 1:N-1
        result[i] = Polyline(l[i], l[i+1])
    end
    result[N] = Polyline(l[N], l[1])
    return result
end

function edges(a::HalfSpace{2,T}) where {T}
    l = Polyline(a.vertex, a.vertex |> translate((a.normal |> rotate(90))...))
    return (Polyline(parameterized(l, -1e6), parameterized(l, 1e6)),)
end

"""
    normals(shape)
"""
function normals(l::Polyline{N,2,T}) where {N,T}
    result = Base.Vector{SVector{2,T}}(undef, N - 1)
    for i in 1:N-1
        result[i] = normalize(vector(l[i], l[i+1]) |> rotate(-90))
    end
    return result
end

function normals(l::Polygon{N,2,T}) where {N,T}
    result = Base.Vector{SVector{2,T}}(undef, N)
    for i in 1:N-1
        result[i] = normalize(vector(l[i], l[i+1]) |> rotate(-90))
    end
    result[N] = normalize(vector(l[N], l[1]) |> rotate(-90))
    return result
end

normals(a::HalfSpace{D,T}) where {D,T} = (a.normal,)

##################################################
# Distances
##################################################

"""
    parameterized(line, a)
"""
parameterized(l::LineSegment{1}, a) = Vertex(
    l[1].x + (l[2].x - l[1].x)*a)
parameterized(l::LineSegment{2}, a) = Vertex(
    l[1].x + (l[2].x - l[1].x)*a,
    l[1].y + (l[2].y - l[1].y)*a)
parameterized(l::LineSegment{3}, a) = Vertex(
    l[1].x + (l[2].x - l[1].x)*a,
    l[1].y + (l[2].y - l[1].y)*a,
    l[1].z + (l[2].z - l[1].z)*a)
parameterized(l::LineSegment{D}, a) where {D} = Vertex(
    l[1][i] + (l[2][i] - l[1][i])*a for i in 1:D)

parameterrange(::LineSegment{D,T}) where {D,T} = (T(0), T(1))

"""
    projection(point, line)

Return the projection of `point` onto `line`.  If the projected point
does not lie on `line`, return `nothing`.

"""
function projection(p::Vertex{D}, l::LineSegment{D}) where {D}
    amin, amax = parameterrange(l)
    d = vector(l)
    a = d ⋅ vector(l[1], p) / (d ⋅ d)
    amin <= a <= amax && return parameterized(l, a)
    return nothing
end

"""
    projection(shape, axis)

Return `shape` projected onto `axis`.

"""
function projection(s::AbstractShape{D}, a::SVector{D}) where {D}
    xmin, xmax = Inf, -Inf
    for v in vertices(s)
        x = a ⋅ vector(v)
        x < xmin && (xmin = x)
        x > xmax && (xmax = x)
    end
    return Polyline(Vertex(xmin), Vertex(xmax))
end

function projection(s::HalfSpace{D,T}, a::SVector{D,T}) where {D,T}
    if iszero(a × s.normal)
        x = a ⋅ vector(s.vertex)
        if a ⋅ s.normal >= 0
            return Polyline(Vertex(T(-Inf)), Vertex(x))
        end
        return Polyline(Vertex(x), Vertex(T(Inf)))
    end
    return Polyline(Vertex(T(-Inf)), Vertex(T(Inf)))
end

"""
    centroid(a, bs...)

Return the centroid (also known as the geometric center) of the
vertices of the shapes, `a, bs...` with equal weight given to each
vertex.

"""
function centroid(a::AbstractShape{1,T}) where {T}
    x = zero(T)
    N = length(vertices(a))
    for v in vertices(a)
        x += v.x
    end
    return Vertex(x/N)
end

function centroid(a::AbstractShape{2,T}) where {T}
    x, y = zero(T), zero(T)
    N = length(vertices(a))
    for v in vertices(a)
        x += v.x
        y += v.y
    end
    return Vertex(x/N, y/N)
end

function centroid(a::AbstractShape{3,T}) where {T}
    x, y, z = zero(T), zero(T), zero(T)
    N = length(vertices(a))
    for v in vertices(a)
        x += v.x
        y += v.y
        z += v.z
    end
    return Vertex(x/N, y/N, z/N)
end

function centroid(a::AbstractShape{D,T}) where {D,T}
    x = fill(zero(T), D)
    N = length(vertices(a))
    for v in vertices(a), i in 1:D
        x[i] += v[i]
    end
    for i in 1:D
        x[i] /= N
    end
    return Vertex(x)
end

function centroid(a::AbstractShape{1}, b::AbstractShape{1})
    aN, bN = length(vertices(a)), length(vertices(b))
    av, bv = centroid(a), centroid(b)
    N = aN + bN
    return Vertex((aN*av.x + bN*bv.x)/N)
end

function centroid(a::AbstractShape{2}, b::AbstractShape{2})
    aN, bN = length(vertices(a)), length(vertices(b))
    av, bv = centroid(a), centroid(b)
    N = aN + bN
    return Vertex((aN*av.x + bN*bv.x)/N, (aN*av.y + bN*bv.y)/N)
end

function centroid(a::AbstractShape{3}, b::AbstractShape{3})
    aN, bN = length(vertices(a)), length(vertices(b))
    av, bv = centroid(a), centroid(b)
    N = aN + bN
    return Vertex((aN*av.x + bN*bv.x)/N, (aN*av.y + bN*bv.y)/N, (aN*av.z + bN*bv.z)/N)
end

function centroid(a::AbstractShape{D,T}, b::AbstractShape{D,T}) where {D,T}
    aN, bN = length(vertices(a)), length(vertices(b))
    av, bv = centroid(a), centroid(b)
    N = aN + bN
    x = fill(zero(T), D)
    for i in 1:D
        x[i] += (aN*av[i] + bN*bv[i])/N
    end
    return Vertex(x)
end

function centroid(shapes::AbstractShape{D,T}...) where {D,T}
    x = fill(zero(T), D)
    N = 0
    for v in flatten(map(vertices, shapes))
        N += 1
        for i in 1:D
            x[i] += v[i]
        end
    end
    for i in 1:D
        x[i] /= N
    end
    return Vertex(x)
end

"""
    smallestradius(shape, vertex)

Return the radius of the smallest circle centered at `vertex` that
contains all vertices in `shape`.

"""
smallestradius(s::AbstractShape{D}, p::Vertex{D}) where {D} =
    sqrt(maximum(v -> sum((v[i] - p[i])^2 for i in 1:D), vertices(s)))

function supportvertex(s::AbstractShape{D,T}, a::SVector{D}) where {D,T}
    xmin = Inf
    p = Vertex(zero(T) for _ in 1:D)
    for v in vertices(s)
        x = a ⋅ vector(v)
        x == xmin && (p = centroid(v, p))
        x < xmin && (xmin = x; p = v)
    end
    return p
end

##################################################
# Transformations
##################################################

struct Transformation{D2,T}
    data::NTuple{D2,T}
end
Transformation(data...) = Transformation(promote(data...))

Base.getindex(a::Transformation, idx) = getindex(Core.getfield(a, :data), idx)

(a::Transformation{4})(b::Transformation{4}) = Transformation(
    a[1]*b[1] - a[2]*b[2],
    a[1]*b[2] + a[2]*b[1],
    a[3] + a[1]*b[3] - a[2]*b[4],
    a[4] + a[2]*b[3] + a[1]*b[4])
(t::Transformation{4})(p::Vertex{2}) = Vertex(t[3] + t[1]*p.x- t[2]*p.y, t[4] + t[1]*p.y + t[2]*p.x)
(t::Transformation{4})(v::SVector{2}) = SVector(t[3] + t[1]*v[1]- t[2]*v[2], t[4] + t[1]*v[2] + t[2]*v[1])
(t::Transformation)(l::Polyline) = Polyline(t(p) for p in l)
(t::Transformation)(l::Polygon) = Polygon(t(p) for p in l)
(t::Transformation)(a::HalfSpace{2}) = HalfSpace(
    t(a.vertex),
    SVector(t[1]*a.normal[1]- t[2]*a.normal[2], t[1]*a.normal[2] + t[2]*a.normal[1]))

"""
    translate(x, y)
    translate(shape, x, y)
"""
translate(x, y) = translate(promote(x, y)...)
translate(x::T, y::T) where T = Transformation(one(T), zero(T), x, y)

"""
    rotate(θ)
    rotate(shape, θ)
"""
rotate(θ) = Transformation(cosd(θ), sind(θ), 0., 0.)

"""
    trace(shape, line)

Return an unsorted vector of parameters, corresponding to each
intersection of `line` with `shape`.

"""
function trace(p::Vertex{2}, l::LineSegment{2})
    b = vector(l)
    c = vector(l[1], p)
    # If `d` is zero, the point is in the line.
    if iszero(b × c)
        # We project `p` onto `l`,
        a = (c ⋅ b) / (b ⋅ b)
        # and check if it is outside the bounds:
        amin, amax = parameterrange(l)
        (a <= amin || amax <= a || isnan(a)) && return ()
        # otherwise, we have a match!
        return (a,)
    end
    return ()
end

function trace(l1::LineSegment{2}, l2::LineSegment{2})
    b = vector(l1)
    c = vector(l2)
    k = vector(l1[1], l2[1])
    d = b × c
    # `t1/d` and `t2/d` are the parameters for the intersection point
    # using `l1` or `l2` respectively.  `t2` is defined a little later.
    t1 = k × c
    # For parallel (or anti–parallel) lines, `d` is zero.
    if iszero(d)
        # If `t1` (or `t2`) is also zero, then the two lines line up.
        if iszero(t1)
            # In this case, we trace the vertices of `l1` by
            # projecting them onto `l2`:
            v1, v2 = vertices(l1)
            n = (c ⋅ c)
            a1 = vector(l2[1], v1) ⋅ c / n
            a2 = vector(l2[1], v2) ⋅ c / n
            # Sort the parameters
            a1, a2 = min(a1, a2), max(a1, a2)
            amin2, amax2 = parameterrange(l2)
            # In the 1D projection, we check if there is no overlap,
            (a2 <= amin2 || amax2 <= a1 || isnan(a1) || isnan(a2)) && return ()
            # in which case we return the parameters corresponding to
            # the interval of  overlap of the line segments:
            return (max(a1, amin2), min(a2, amax2))
        end
        # If the lines do not line up, then they do not intersect.
        return ()
    end
    # For non–parallel lines, we check if the parameter `t1/d` is in
    # `l1`, and return `t2/d` if it is in `l2` (TODO: include simple
    # proof why this works).
    t2 = k × b
    amin1, amax1 = parameterrange(l1)
    amin2, amax2 = parameterrange(l2)
    if amin1 <= t1/d <= amax1 && amin2 <= t2/d <= amax2
        return (t2/d,)
    end
    return ()
end

trace(s::AbstractLineShape, l::LineSegment) = flatten(map(e -> trace(e, l), edges(s)))

"""
    intersects(shape1, shape2)
"""
intersects(p1::Vertex, p2::Vertex) = p1 == p2

# See `trace` above for detailed comments of the code.
function intersects(p::Vertex{2}, l::LineSegment{2})
    b = vector(l)
    c = vector(l[1], p)
    if iszero(b × c)
        a = (c ⋅ b) / (b ⋅ b)
        amin, amax = parameterrange(l)
        (a <= amin || amax <= a || isnan(a)) && return false
        return true
    end
    return false
end

# See `trace` above for detailed comments of the code.
function intersects(l1::LineSegment{2}, l2::LineSegment{2})
    b = vector(l1)
    c = vector(l2)
    k = vector(l1[1], l2[1])
    d = b × c
    t1 = k × c
    if iszero(d)
        if iszero(t1)
            v1, v2 = vertices(l1)
            n = (c ⋅ c)
            a1 = vector(l2[1], v1) ⋅ c / n
            a2 = vector(l2[1], v2) ⋅ c / n
            a1, a2 = min(a1, a2), max(a1, a2)
            amin2, amax2 = parameterrange(l2)
            (a2 <= amin2 || amax2 <= a1 || isnan(a1) || isnan(a2)) && return false
            return true
        end
        return false
    end
    t2 = k × b
    amin1, amax1 = parameterrange(l1)
    amin2, amax2 = parameterrange(l2)
    if amin1 <= t1/d <= amax1 && amin2 <= t2/d <= amax2
        return true
    end
    return false
end

intersects(a::AbstractLineShape{N,D}, b::Vertex{D}) where {N,D} = any(e->intersects(e, b), edges(a))
intersects(a::AbstractLineShape{N,D}, b::LineSegment{D}) where {N,D} = any(e->intersects(e, b), edges(a))
intersects(a::Union{Vertex,LineSegment}, b::AbstractLineShape) = intersects(b, a)

const product = Iterators.product
function intersects(a::AbstractLineShape{N1,D}, b::AbstractLineShape{N2,D}) where {N1,N2,D}
    !intersects(BoundingBox(a), BoundingBox(b)) && return false
    for ae in edges(a), be in edges(b)
        intersects(ae, be) && return true
    end
    return false
end

"""
    inside(shape1, shape2, θ)

Determine if `shape2` is inside `shape1`, using a rays from each point
at an angle `θ`.  A shape is determined to be inside if either:

 * the two shapes intersect, or;

 * a ray starting at one of the vertices of `shape2` intersects
   `shape1` an odd number of times.

"""
function inside(l::AbstractLineShape{N,2}, p::Vertex{2}) where N
    # Ray from p through shape center
    c = centroid(l)
    p2c = normalize(c .- p)
    ray_endpoint = Vertex((p .+ 1e6 * p2c)...)
    return inside(l, p, ray_endpoint)
end
function inside(l::AbstractLineShape{N,2}, p::Vertex{2}, θ::Number) where N
    ray_endpoint = p |> translate(cosd(θ)*1e6, sind(θ)*1e6)
    return inside(l, p, ray_endpoint)
end
function inside(l::AbstractLineShape{N,2}, p1::Vertex{2}, p2::Vertex{2}) where N
    ray = Polyline(p1, p2)
    n = count(edges(l)) do e
        if intersects(e, ray)
            i1 = intersects(e[1], ray)
            i2 = intersects(e[2], ray)
            # Intersecting both end points counts as a double
            # intersection.
            i1 && i2 && return false
            # Not intersecting either end point is a valid
            # intersection.
            !(i1 || i2) && return true
            # Intersecting one vertex only counts if the the other
            # vertex of the edge lies clockwise w.r.t. to the ray.
            # This overcomes the issue of tracing corners
            d = vector(ray)
            (i1 && d × vector(ray[1], e[2]) >= 0) && return true
            (i2 && d × vector(ray[1], e[1]) >= 0) && return true
        end
        return false
    end
    return isodd(n)
end

"""
    separatingaxis(a, b)

Return a tuple of `(d, axis)`, where `axis` satisfies the condition
that the projections of `a` and `b` do not overlap with distance `-d`
between them.  If no such axis is found, return instead the `axis`
such that `d*axis` is the minimum translation vector applied to `b`
required to push `b` out of `a`.

"""
function separatingaxis(a::AbstractShape{D,T}, b::AbstractShape{D,T}) where {D,T}
    d = Inf
    axis = zero(SVector{D,T})
    for n in normals(a)
        an = projection(a, n)
        bn = projection(b, n)
        d′ = an[2].x - bn[1].x
        if d′ < d
            d = d′
            axis = n
        end
        d′ = bn[2].x - an[1].x
        if d′ < d
            d = d′
            axis = -n
        end
        d < 0 && return (d, axis)
    end
    return d, axis
end

end #module


text:
### Generic MUSE code

abstract type AbstractMuseProblem end

## interface to be implemented by specific problem types

function ∇θ_logLike end
function logLike_and_∇z_logLike end
function sample_x_z end
logPriorθ(prob::AbstractMuseProblem, θ) = 0
standardizeθ(prob::AbstractMuseProblem, θ) = θ

# this can also be overriden by specific problems
# the default does LBFGS using the provided logLike_and_∇z_logLike
function ẑ_at_θ(prob::AbstractMuseProblem, x, z₀, θ; ∇z_logLike_atol)
    soln = optimize(Optim.only_fg(z -> .-logLike_and_∇z_logLike(prob, x, z, θ)), z₀, Optim.LBFGS(), Optim.Options(g_tol=∇z_logLike_atol))
    soln.minimizer, soln
end

### MUSE result

@doc doc"""

Stores the result of a MUSE run. Can be constructed by-hand as
`MuseResult()` and passed to any of the inplace `muse!`, `get_J!`, or
`get_H!`.

Fields:

* `θ` — The estimate of the $\theta$ parameters.
* `Σ, Σ⁻¹` — The approximate covariance of $\theta$ and its inverse.
* `H, J` — The $H$ and $J$ matrices which form the covariance (see
  [Millea & Seljak, 2021](https://arxiv.org/abs/2112.09354))
* `gs` — The MAP gradient sims used to compute `J`.
* `Hs` — The jacobian sims used to compute `H`.
* `dist` — A `Normal` or `MvNormal` built from `θ` and `Σ`, for
  convenience.
* `history` — Internal diagnostic info from the run.
* `rng` — RNG used to generate sims for this run (so the same sims can
  be reused if resuming later).
* `time` — Total `Millisecond` wall-time spent computing the result.

"""
Base.@kwdef mutable struct MuseResult
    θ = nothing
    H = nothing
    J = nothing
    Σ⁻¹ = nothing
    Σ = nothing
    dist = nothing
    history = []
    gs = []
    Hs = []
    rng = nothing
    time = Millisecond(0)
end

### MUSE solver

@doc doc"""

    muse(prob::AbstractMuseProblem, θ₀; kwargs...)
    muse!(result::MuseResult, prob::AbstractMuseProblem, [θ₀=nothing]; kwargs...)

Run the MUSE estimate. The `muse!` form resumes an existing result. If the
`muse` form is used instead, `θ₀` must give a starting guess for $\theta$.

See [`MuseResult`](@ref) for description of return value.

Optional keyword arguments:

* `rng` — Random number generator to use. Taken from `result.rng` or `Random.default_rng()` if not passed.
* `z₀` — Starting guess for the latent space MAP.
* `maxsteps = 50` — Maximum number of iterations.
* `θ_rtol = 1e-1` — Error tolerance on $\theta$ relative to its standard deviation.
* `∇z_logLike_atol = 1e-2` — Absolute tolerance on the $z$-gradient at the MAP solution.
* `nsims = 100` — Number of simulations.
* `α = 0.7` — Step size for root-finder.
* `progress = false` — Show progress bar.
* `pmap` — Parallel map function.
* `regularize = identity` — Apply some regularization after each step.
* `H⁻¹_like = nothing` — Initial guess for the inverse Jacobian of $s^{\rm MUSE}(\theta)$
* `H⁻¹_update` — How to update `H⁻¹_like`. Should be `:sims`, `:broyden`, or `:diagonal_broyden`.
* `broyden_memory = Inf` — How many past steps to keep for Broyden updates.
* `checkpoint_filename = nothing` — Save result to a file after each iteration.
* `get_covariance = false` — Also call `get_H` and `get_J` to get the full covariance.

"""
muse(args...; kwargs...) = muse!(MuseResult(), args...; kwargs...)

function muse!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = nothing;
    rng = nothing,
    z₀ = nothing,
    maxsteps = 50,
    θ_rtol = 1e-1,
    ∇z_logLike_atol = 1e-2,
    nsims = 100,
    α = 0.7,
    progress = false,
    pmap = _map,
    batch_size = 1,
    regularize = identity,
    H⁻¹_like = nothing,
    H⁻¹_update = :sims,
    broyden_memory = Inf,
    checkpoint_filename = nothing,
    get_covariance = false
)

    rng = @something(rng, result.rng, copy(Random.default_rng()))
    θunreg = θ = θ₀ = standardizeθ(prob, @something(result.θ, θ₀))
    z₀ = @something(z₀, sample_x_z(prob, copy(rng), θ₀).z)
    local H⁻¹_post, g_like_sims
    history = result.history

    result.rng = _rng = copy(rng)
    xz_sims = [sample_x_z(prob, _rng, θ) for i=1:nsims]
    xs = [[prob.x];  getindex.(xz_sims, :x)]
    ẑs = [[z₀];      getindex.(xz_sims, :z)]

    # set up progress bar
    pbar = progress ? RemoteProgress((maxsteps-length(result.history))*(nsims+1)÷batch_size, 0.1, "MUSE: ") : nothing

    try

        for i = (length(result.history)+1):maxsteps

            t₀ = now()

            if i > 1
                _rng = copy(rng)
                xs = [[prob.x]; [sample_x_z(prob, _rng, θ).x for i=1:nsims]]
            end

            if i > 2
                Δθ = history[end].θ - history[end-1].θ
                norm(Δθ ./ θ) < θ_rtol && break
            end

            # MUSE gradient
            gẑs = pmap(xs, ẑs, fill(θ,length(xs)); batch_size) do x, ẑ_prev, θ
                local ẑ, history = ẑ_at_θ(prob, x, ẑ_prev, θ; ∇z_logLike_atol)
                g = ∇θ_logLike(prob, x, ẑ, θ)
                progress && ProgressMeter.next!(pbar)
                (;g, ẑ, history)
            end
            ẑs = getindex.(gẑs, :ẑ)
            ẑ_history_dat, ẑ_history_sims = peel(getindex.(gẑs, :history))
            g_like_dat, g_like_sims = peel(getindex.(gẑs, :g))
            g_like = g_like_dat .- mean(g_like_sims)
            g_prior = AD.gradient(AD.ForwardDiffBackend(), θ -> logPriorθ(prob, θ), θ)[1]
            g_post = g_like .+ g_prior

            # Jacobian
            h⁻¹_like_sims = -1 ./ var(collect(g_like_sims))
            H⁻¹_like_sims = h⁻¹_like_sims isa Number ? h⁻¹_like_sims : Diagonal(h⁻¹_like_sims)
            if (H⁻¹_like == nothing) || (H⁻¹_update == :sims)
                H⁻¹_like = H⁻¹_like_sims
            elseif i > 2 && (H⁻¹_update in [:broyden, :diagonal_broyden])
                # on subsequent steps, do a Broyden's update using at
                # most the previous `broyden_memory` steps
                j₀ = Int(max(2, i - broyden_memory))
                H⁻¹_like = history[j₀-1].H⁻¹_like_sims
                for j = j₀:i-1
                    Δθ      = history[j].θ      - history[j-1].θ
                    Δg_like = history[j].g_like - history[j-1].g_like
                    H⁻¹_like = H⁻¹_like + ((Δθ - H⁻¹_like * Δg_like) / (Δθ' * H⁻¹_like * Δg_like)) * Δθ' * H⁻¹_like
                    if H⁻¹_update == :diagonal_broyden
                        H⁻¹_like = Diagonal(H⁻¹_like)
                    end
                end
            end

            H_prior = AD.hessian(AD.ForwardDiffBackend(), θ -> logPriorθ(prob, θ), θ)[1]
            H⁻¹_post = inv(inv(H⁻¹_like) + H_prior)

            t = now() - t₀
            push!(
                history,
                (;θ, θunreg,
                g_like_dat, g_like_sims, g_like, g_prior, g_post,
                H⁻¹_post, H_prior, H⁻¹_like, H⁻¹_like_sims,
                ẑ_history_dat, ẑ_history_sims, t)
            )

            # Newton-Rhapson step
            θunreg = θ .- α .* (H⁻¹_post * g_post)
            θ = regularize(θunreg)

            (checkpoint_filename != nothing) && save(checkpoint_filename, "result", result)

        end

    finally

        progress && ProgressMeter.finish!(pbar)

    end

    result.time += sum(getindex.(history,:t))
    result.θ = θunreg
    result.gs = collect(g_like_sims)
    if get_covariance
        get_J!(result, prob)
        get_H!(result, prob)
    end
    result

end

function get_H!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = result.θ;
    fdm :: FiniteDifferenceMethod = central_fdm(3,1),
    ∇z_logLike_atol = 1e-8,
    rng = Random.default_rng(),
    nsims = 10,
    step = nothing,
    pmap = _map,
    batch_size = 1,
    pmap_over = :auto,
    progress = false,
    skip_errors = false,
)

    θ₀ = standardizeθ(prob, @something(θ₀, result.θ))
    nsims_remaining = nsims - length(result.Hs)
    (nsims_remaining <= 0) && return
    pbar = progress ? RemoteProgress(nsims_remaining*(1+length(θ₀))÷batch_size, 0.1, "get_H: ") : nothing
    t₀ = now()

    # generate simulation locally, advancing rng, and saving rng state to be reused remotely
    xs_zs_rngs = map(1:nsims_remaining) do i
        _rng = copy(rng)
        (x, z) = sample_x_z(prob, rng, θ₀)
        (x, z, _rng)
    end

    # initial fit at fiducial, used at starting points for finite difference below
    ẑ₀s_rngs = pmap(xs_zs_rngs; batch_size) do (x, z, rng)
        ẑ, = ẑ_at_θ(prob, x, z, θ₀; ∇z_logLike_atol)
        progress && ProgressMeter.next!(pbar)
        (ẑ, rng)
    end

    # finite difference Jacobian
    pmap_sims, pmap_jac = (pmap_over == :jac || (pmap_over == :auto && length(θ₀) > nsims_remaining)) ? (_map, pmap) : (pmap, _map)
    append!(result.Hs, skipmissing(pmap_sims(ẑ₀s_rngs; batch_size) do (ẑ₀, rng)
        try
            return first(pjacobian(fdm, θ₀, step; pmap=pmap_jac, batch_size, pbar) do θ
                x, = sample_x_z(prob, copy(rng), θ)
                ẑ, = ẑ_at_θ(prob, x, ẑ₀, θ₀; ∇z_logLike_atol)
                ∇θ_logLike(prob, x, ẑ, θ₀)
            end)
        catch err
            if skip_errors && !(err isa InterruptException)
                @warn err
                return missing
            else
                rethrow(err)
            end
        end
    end))

    result.H = (θ₀ isa Number) ? mean(first.(result.Hs)) :  mean(result.Hs)
    result.time += now() - t₀
    finalize_result!(result, prob)

end

function get_J!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = nothing;
    ∇z_logLike_atol = 1e-1,
    rng = Random.default_rng(),
    nsims = 100,
    pmap = _map,
    batch_size = 1,
    progress = false,
    skip_errors = false,
    covariance_method = LinearShrinkage(target=DiagonalCommonVariance(), shrinkage=:rblw),
)

    θ₀ = standardizeθ(prob, @something(θ₀, result.θ))
    nsims_remaining = nsims - length(result.gs)

    if nsims_remaining > 0

        pbar = progress ? RemoteProgress(nsims_remaining÷batch_size, 0.1, "get_J: ") : nothing

        (xs, zs) = map(Base.vect, map(1:nsims_remaining) do i
            sample_x_z(prob, rng, θ₀)
        end...)

        append!(result.gs, skipmissing(pmap(xs, zs, fill(θ₀,length(xs)); batch_size) do x, z, θ₀
            try
                ẑ, = ẑ_at_θ(prob, x, z, θ₀; ∇z_logLike_atol)
                g = ∇θ_logLike(prob, x, ẑ, θ₀)
                progress && ProgressMeter.next!(pbar)
                return g
            catch err
                if skip_errors && !(err isa InterruptException)
                    @warn err
                    return missing
                else
                    rethrow(err)
                end
            end
        end))

    end

    result.J = (θ₀ isa Number) ? var(result.gs) : cov(covariance_method, identity.(result.gs))
    finalize_result!(result, prob)

end

function finalize_result!(result::MuseResult, prob::AbstractMuseProblem)
    @unpack H, J = result
    if H != nothing && J != nothing
        H_prior = -AD.hessian(AD.ForwardDiffBackend(), θ -> logPriorθ(prob, θ), result.θ)[1]
        result.Σ⁻¹ = H' * inv(J) * H + H_prior
        result.Σ = inv(result.Σ⁻¹)
        if length(result.θ) == 1
            result.dist = Normal(result.θ[1], sqrt(result.Σ[1]))
        else
            result.dist = MvNormal(result.θ, result.Σ)
        end
    end
    result
end

ext: jl


text:
using GNSSTools
using FFTW
using AllanDeviations
using Statistics
# using LsqFit
# using NNLS
using PyPlot
pygui(true)

# function psd_at_f(h_parms, f, v_0)
#     h₋₂, h₋₁, h₀, h₁, h₂ = h_parms
#     return v_0^2 * (h₋₂*f^(-4) + h₋₁*f^(-3) + h₀*f^(-2) + h₁*f^(-1) + h₂*f^0)
# end

# function fit_h_parms_to_signal_PSD(psd_dBc, f_c=10e6,
#                                    freqs=[1,10,100,1e3,1e4,1e5])
#     N = length(freqs)
#     @assert N == length(psd_dBc)
#     b = zeros(N-1)
#     A = zeros(N-1)
#     for i in 1:N-1
#         j = N - i
#         A[i] = (psd_dBc[j] + psd_dBc[j+1])/2 + 10*log10(freqs[j+1] - freqs[j])
#         if i == 1
#             b[i] = 2*10^(A[i]/10) / (freqs[j+1] - freqs[j])  # b₀
#         elseif i == 2
#             b[i] = 2*10^(A[i]/10) / log(freqs[j+1] - freqs[j])  # b₋₁
#         else
#             b[i] = 2*10^(A[i]/10) / ((2-i)*(freqs[j+1]^(2-i) - freqs[j]^(2-i)))
#         end
#     end
#     # order is h2, h1, h0, h-1, h-2
#     h = b ./ f_c.^2
#     # reverse order to h-2, h-1, h₀, h₁, h₂
#     return reverse(h)
# end

function calc_Sy(x, f_s)
    N = length(x)
    X = fft(x)[1:floor(Int, N/2)+1]
    X = (1/(f_s*N)) .* abs2.(X)
    X[2:end-1] .= 2 .* X[2:end-1]
    X_log10 = 10 .* log10.(X)
    frequencies = Array(range(0, f_s/2, length=length(X)))
    return (X, X_log10, frequencies)
end

# f_s = 36e6
f_carrier1 = 16.3676e6
f_carrier2 = 10e6
f_s = 1*f_carrier1
t_length = 1
N = floor(Int, t_length*f_s)
N_over_2 = floor(Int, N/2)
t = Array(range(0, t_length, length=N))
frequencies = Array(range(0, f_s/2, length=N_over_2+1))
# carrier1 = cis.(2π.*f_carrier1.*t)
# carrier2 = cis.(2π.*f_carrier2.*t)
phase_noise = zeros(Complex{Float64}, N)
# signal = zeros(Complex{Float64}, N)
dBc_freq = [1., 1e1, 1e2, 1e3, 1e4]
X = zeros(Complex{Float64}, N_over_2+1)
x = zeros(Complex{Float64}, N)
rakon_it5300B = Dict("h_parms"=>[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15],
                     "dBc"=>[-57, -88, -112, -130, -140],
                     "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                     "name"=>"TCXO, Rakon IT5300B",
                     "f_0"=>f_carrier1)
agilent_e4424b = Dict("h_parms"=>[2.8e-9, 1.8e-11, 1.2e-10, 1e-11, 1.8e-15],
                      "dBc"=>[-88, -118, -128, -135, -147],
                      "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                      "name"=>"OCXO, Agilent E4424B",
                      "f_0"=>10e6)
isotemp_91_1 = Dict("h_parms"=>[2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16],
                    "dBc"=>[-89, -120, -140, -151, -154],
                    "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                    "name"=>"OCXO, ISOTEMP 91-1",
                    "f_0"=>f_carrier2)
oscilloquartz_8607 = Dict("h_parms"=>[3.9e-14, 8.7e-13, 2.2e-14, 3.2e-13, 6.4e-15],
                          "dBc"=>[-122, -137, -143, -145, -145],
                          "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                          "name"=>"OCXO, Oscilloquartz 8607",
                          "f_0"=>f_carrier2)
srs_prs10 = Dict("h_parms"=>[5.3e-11, 4.7e-11, 1.6e-15, 9.6e-14, 9.9e-16],
                 "dBc"=>[-103, -135, -150, -152, -153],
                 "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                 "name"=>"Rubidum, SRS PRS10",
                 "f_0"=>f_carrier2)
cesium = Dict("h_parms"=>[1.5e-12, 6.2e-11, 7.6e-23, 1.4e-14, 6.2e-19],
              "dBc"=>[-105, -135, -160, -170, -190],
              "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
              "name"=>"Cesium, Typical Values",
              "f_0"=>f_carrier2)
mti_250l = Dict("h_parms"=>[1.7e-9, 1e-10, 1.6e-10, 3.5e-13, 1.8e-16],
              "dBc"=>[-90, -120, -140, -155, -160],
              "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
              "name"=>"OCXO, MTI 250L",#, Low g-sensitivity",
              "f_0"=>f_carrier2)

oscillators = [rakon_it5300B, agilent_e4424b, isotemp_91_1, oscilloquartz_8607,
               srs_prs10, cesium, mti_250l]
# figsize = (6.5, 4)
colorz = ["k", "b", "g", "r", "y", "m", "c", "gray"]
fig = figure(figsize=(6.5,4))
ax1 = fig.add_subplot(1,1,1)
legend_markers = []
marker_names = []
for i in 1:length(oscillators)
# for i in 1:1
    color = colorz[i]
    oscillator = oscillators[i]
    if oscillator["name"] == "TCXO, Rakon IT5300B"
        # carrier = carrier1
        f_carrier = f_carrier1
    else
        # carrier = carrier2
        f_carrier = f_carrier2
    end
    h_parms = oscillator["h_parms"] ./ oscillator["f_0"]^2
    generate_phase_noise!(phase_noise, t_length, f_carrier, h_parms)
    # ϕ = real.(phase_noise)
    # signal = carrier .* cis.(ϕ)
    # signal .= carrier .* cis.(real.(phase_noise))
    # x = signal
    x .= real.(phase_noise) .+ 0im
    # X, X_log, freqs = calc_Sy(x, f_s)
    fft!(x);
    X .= abs2.(x[1:N_over_2+1]) / (f_s*N);
    X[2:end-1] .= 2 .* X[2:end-1];
    X .= 10 .* log10.(X);
    # X_max = maximum(real.(X));
    # X = X .- X_max;
    # ax1.plot(frequencies[2:end] .- f_carrier, X[2:end], color,
    #          label=string(oscillator["name"], " from h Parms"))
    ax1.plot(frequencies[2:end], real.(X[2:end]), color)#),
            #  label=oscillator["name"])
    # ax1.plot(frequencies[2:end], real.(X[2:end]), color,
            #  label=string(oscillator["name"], " from h Parms"))
    ax1.scatter(dBc_freq, oscillator["dBc"].+(10*log10(2)),
                c=color)#, label=string(oscillator["name"], " Truth"))
    legend_marker = plt.Line2D((0,1),(0,0), color=color, marker="o", linestyle="-")
    push!(legend_markers, legend_marker)
    push!(marker_names, oscillator["name"])
end
xscale("log")
xlim([1, 1e6])
xlabel("Frequency (Hz)")
ylabel("dBc/Hz")
subplots_adjust(top=0.96, left=0.12, right=0.88, bottom=0.13)
ax1.legend(legend_markers, marker_names)

# savefig("figures/ch2_oscillator_psds.pdf", dpi=300)


text:
#= test/trapezoidal.jl
=#

@testset "trapezoidal" begin
  m = 20
  n = 10
  data_up = rand(ComplexF64, m, n)
  data_lo = data_up'

  @time for (t, data, trilu) in ((LowerTrapezoidal, data_up, tril),
                                 (UpperTrapezoidal, data_lo, triu))
    for T in (Float32, Float64, ComplexF32, ComplexF64)
      let n, m
          println("  $t/$T")

          A = t{T}(T <: Real ? real(data) : data)
          m, n = size(A)
          F = Matrix(A)
          @test trilu(A.data) ≈ F

          xm = rand(T, m)
          xn = rand(T, n)
          @test A  *xn ≈ F  *xn
          @test A' *xm ≈ F' *xm
          @test transpose(A)*xm ≈ transpose(F)*xm

          Bm = rand(T, m, m)
          Bn = rand(T, n, n)
          @test A   *Bn   ≈ F   *Bn
          @test A   *Bn'  ≈ F   *Bn'
          @test A   *transpose(Bn) ≈ F   *transpose(Bn)
          @test A'  *Bm   ≈ F'  *Bm
          @test A'  *Bm'  ≈ F'  *Bm'
          @test transpose(A) *Bm   ≈ transpose(F) *Bm
          @test transpose(A) *transpose(Bm) ≈ transpose(F) *transpose(Bm)
          @test Bm  *A    ≈ Bm  *F
          @test Bn  *A'   ≈ Bn  *F'
          @test Bn  *transpose(A)  ≈ Bn  *transpose(F)
          @test Bm' *A    ≈ Bm' *F
          @test Bn' *A'   ≈ Bn' *F'
          @test transpose(Bm)*A    ≈ transpose(Bm)*F
          @test transpose(Bn)*transpose(A)  ≈ transpose(Bn)*transpose(F)
      end
    end
  end
end


text:
using Jui
using Test

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


text:
function is_sense_legit(testSense::Symbol)::Bool

    # we allow buy and sell operations -
    approvedSet = Set{Symbol}()
    push!(approvedSet,:buy)
    push!(approvedSet,:sell)

    # check - is the testSense a member of the approvedSet?
    if (in(testSense,approvedSet) == true)
        return true
    end

    return false
end

function is_positive_value(value::T where T<:Number)::Bool

    # we should have only positive values -
    if (sign(value)==1.0 || sign(value) == 1 || sign(value) == 0 || sign(value) == 0.0)
        return true
    end

    return false
end


text:
const Terminal = Union{Symbol, Function}

"""
Implementation of Koza-type (tree-based) Genetic Programming

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `initialization`: A strategy for population initialization (default: `:half`)
  - Possible values: `:grow`, `:full`, `:half`

"""
@kwdef struct TreeGP <: AbstractOptimizer
    populationSize::Integer = 50
    terminals::Dict{Terminal, Int} = Dict(:x=>1, rand=>1)
    functions::Dict{Function, Int} = Dict( f=>2 for f in [+,-,*,/] )
    mindepth::Int = 0
    maxdepth::Int = 10
    initialization::Symbol = :grow
    optimizer::AbstractOptimizer = GA()
end
function TreeGP(pop::Integer, term::Vector{Terminal}, func::Vector{Function}; kwargs...)
    terminals = Dict(t=>1 for t in term)
    functions = Dict(f=>2 for f in func)
    TreeGP(;populationSize=pop, terminals=terminals, functions=functions, kwargs...)
end

population_size(method::TreeGP) = method.populationSize
default_options(method::TreeGP) = (iterations=1000, abstol=1e-15)
terminals(m::TreeGP) = Symbol[t for t in keys(m.terminals) if isa(t,Symbol)] |> sort!
function summary(m::TreeGP)
    par = join(terminals(m),",")
    "TreeGP[P=$(m.populationSize),Parameter[$(par)],$(keys(m.functions))]"
end

function rand(t::TreeGP; maxdepth=2, mindepth=maxdepth-1, method=:grow)
    # println("depth: $maxdepth => $mindepth")
    @assert maxdepth > mindepth "`maxdepth` must be larger then `mindepth`"
    tl = length(t.terminals)
    fl = length(t.functions)
    root = if (maxdepth == 0  || ( method == :grow && rand() < tl/(tl+fl) ) ) && mindepth <= 0
        term = rand(keys(t.terminals))
        if isa(term, Symbol)
            term
        else
            dim = t.terminals[term]
            dim == 1 ? rand() : rand(dim)
        end
    else
        rand(keys(t.functions))
    end
    if isa(root, Function)
        args = Any[]
        for i in 1:t.functions[root]
            arg = rand(t, maxdepth=maxdepth-1, mindepth=mindepth-1, method=method)
            # println("arg_$i: $arg ")
            push!(args, arg)
        end
        Expr(:call, root, args...)
    else
        return root
    end
end

"""
    initial_population(method::TreeGP, expr::{Expr,Nothing}=nothing)

Initialize a random population of expressions derived from `expr`.
"""
function initial_population(method::TreeGP, expr::Union{Expr,Nothing}=nothing)
    n = population_size(method)
    return [rand(method, mindepth=method.mindepth, maxdepth=method.maxdepth) for i in 1:n]
end

height(ex) = isa(ex, Expr) ? maximum( height(e) for e in ex.args )+1 : 0
nodes(ex) = !isa(ex, Expr) ? 1 : length(ex.args) > 0 ? sum( nodes(e) for e in ex.args ) : 0
length(ex) = nodes(ex)

function depth(root, ex; d=0)
    return if root == ex
        d
    elseif isa(root, Expr)
        maximum( depth(e, ex, d=d+1) for e in root.args )
    else
        -1
    end
end

function copyto!(ex1::Expr, ex2::Expr)
    ex1.head = ex2.head
    ex1.args = deepcopy(ex2.args)
end

function randsubtree(ex)
    !isa(ex, Expr) && return ex
    csize = map(nodes, ex.args[2:end])
    cidx = rand(1:sum(csize)+1)-1
    if cidx == 0
        ex
    else
        nidx = findfirst(i->cidx<=i, cumsum(csize))
        randsubtree(ex.args[nidx+1])
    end
end

function getindex(ex::Expr, idx::Int)
    if idx == 0
        ex
    else
        csize = cumsum(map(nodes, ex.args[2:end]))
        nidx = findfirst(i->i>=idx, csize)
        nex = ex.args[nidx+1]
        !isa(nex, Expr) ? nex : nex[csize[nidx]-idx]
    end
end

function setindex!(ex::Expr, subex, idx::Int)
    csize = cumsum(map(nodes, ex.args[2:end]))
    nidx = findfirst(i->i>=idx, csize)
    nex = ex.args[nidx+1]
    if !isa(nex, Expr) || csize[nidx]-idx == 0
        ex.args[nidx+1] = subex
    else
        nex[csize[nidx]-idx] = subex
    end
end

function crosstree(t1::Expr, t2::Expr)
    tt1, tt2 = copy(t1), copy(t2)
    i, j = rand(1:nodes(t1)-1), rand(1:nodes(t2)-1)
    ex1 = tt1[i]
    ex2 = tt2[j]
    tt1[i] = ex2
    tt2[j] = ex1
    tt1, tt2
end

function mutatetree(method::TreeGP)
    function mutation(recombinant::Expr)
        i = rand(1:nodes(recombinant)-1)
        th = depth(recombinant, recombinant[i])
        recombinant[i] = rand(method, maxdepth=max(0, method.maxdepth-th))
        recombinant
    end
    return mutation
end

"""
- `fitness`: A collection of fitness cases required for evaluation of a fitness function
"""
mutable struct GPState{T,IT} <: Evolutionary.AbstractOptimizerState
    ga::GAState{T,IT}
end
value(s::GPState) = s.ga.fitness
minimizer(s::GPState) = s.ga.fittest

"""Initialization of GA algorithm state"""
function initial_state(method::TreeGP, options, objfun, population)
    # setup initial state
    return GPState(initial_state(method.optimizer, options, objfun, population))
end

function update_state!(objfun, constraints, state::GPState, population::AbstractVector{IT}, method, itr) where {IT}
    # perform GA step
    res = update_state!(objfun, constraints, state.ga, population, method.optimizer, itr)
    # simplify expressions
    for i in 1:length(population)
        simplify!(population[i])
    end
    return res
end

function symbols(ex::Expr)
    syms = Symbol[]
    for e in ex.args
        isa(e, Symbol) && push!(syms, e)
        isa(e, Expr) && append!(syms, symbols(e))
    end
    unique!(syms)
end

struct Expression
    expr::Expr
    syms::Vector{Symbol}
end
Expression(ex::Expr) = Expression(ex, sort!(symbols(ex)))
show(io::IO, e::Expression) = infix(io, e.expr)
(e::Expression)(val) = evaluate(val, e.expr, e.syms)

function evaluate(val, ex::Expr, psyms::Vector{Symbol})
    exprm = ex.args
    exvals = (isa(nex, Expr) || isa(nex, Symbol) ? evaluate(val, nex, psyms) : nex for nex in exprm[2:end])
    exprm[1](exvals...)
end

function evaluate(val, ex::Symbol, psyms::Vector{Symbol})
    pidx = findfirst(isequal(ex), psyms)
    val[pidx]
end

iszeronum(root) = isa(root, Number) && iszero(root)
isonenum(root) = isa(root, Number) && isone(root)

function simplify!(root)
    # return non-expression
    !isa(root, Expr) && return root
    if root.head == :call
        for i in 2:3
            !isa(root.args[i], Expr) && continue
            root.args[i] = simplify!(root.args[i])
        end
        # simplification rules
        if root.args[1] == (-) && root.args[2] == root.args[3] # additive inverse
            root = 0
        elseif root.args[1] == (/) && (root.args[2] == root.args[3] || iszeronum(root.args[3])) # look for multiplicative inverse
            root = 1
        elseif (root.args[1] == (*) || root.args[1] == (/)) &&  (iszeronum(root.args[2]) || iszeronum(root.args[3])) # look for 0
            root = 0
        elseif (root.args[1] == (+) || root.args[1] == (-)) && iszeronum(root.args[3]) # look for 0
            root = root.args[2]
        elseif root.args[1] == (+) && iszeronum(root.args[2]) # look for 0
            root = root.args[3]
        elseif (root.args[1] == (*) || root.args[1] == (/)) && isonenum(root.args[3]) # x*1 = x || x/1 = x
            root = root.args[2]
        elseif root.args[1] == (*) && isonenum(root.args[2]) # 1*x = x
            root = root.args[3]
        else
            if root.args[1] == (+) && root.args[2] == root.args[3] # + into *
                root.args[1] = (*)
                root.args[2] = 2
            end
            if root.args[1] == (*) && # (z/d)*(x*y) = x*z (if y = d) || y*z (if x = d)
                isa(root.args[2], Expr) && root.args[2].args[1] == (/) &&
                isa(root.args[3], Expr) && root.args[3].args[1] == (*)
                x = root.args[3].args[2]
                y = root.args[3].args[3]
                z = root.args[2].args[2]
                d = root.args[2].args[3]
                if x == d
                    root.args[2] = z
                    root.args[3] = y
                elseif y == d
                    root.args[2] = z
                    root.args[3] = x
                elseif isa(d, Number) && isa(x, Number) # (z/d)*(x*y) = (y*z)*(x/d) = (y*z)*w
                    root.args[3] = x/d
                    root.args[2].args[1] = *
                    root.args[2].args[2] = y
                    root.args[2].args[3] = z
                elseif isa(d, Number) && isa(y, Number)  # (z/d)*(x*y) = (x*z)*(y/d) = (x*z)*w
                    root.args[3] = y/d
                    root.args[2].args[1] = *
                    root.args[2].args[2] = x
                    root.args[2].args[3] = z
                end
            end
            if root.args[1] == (*) && # (x*y)*(z/d) = x*z (if y = d) || y*z (if x = d)
                isa(root.args[2], Expr) && root.args[2].args[1] == (*) &&
                isa(root.args[3], Expr) && root.args[3].args[1] == (/)
                x = root.args[2].args[2]
                y = root.args[2].args[3]
                z = root.args[3].args[2]
                d = root.args[3].args[3]
                if x == d
                    root.args[2] = y
                    root.args[3] = z
                elseif y == d
                    root.args[2] = x
                    root.args[3] = z
                elseif isa(d, Number) && isa(x, Number) # (x*y)*(z/d) = (x/d)*(y*z) = w*(y*z)
                    root.args[2] = x/d
                    root.args[3].args[1] = *
                    root.args[3].args[2] = y
                    root.args[3].args[3] = z
                elseif isa(d, Number) && isa(y, Number)  # (x*y)*(z/d) = (y/d)*(x*z) = w*(x*z)
                    root.args[2] = y/d
                    root.args[3].args[1] = *
                    root.args[3].args[2] = x
                    root.args[3].args[3] = z
                end
            end
            if (root.args[1] == (+) || root.args[1] == (-)) && # x ± (y ± z) = (x ± y) ± z = w ± z
                isa(root.args[2], Number) &&
                isa(root.args[3], Expr) && (root.args[3].args[1] == (+) || root.args[3].args[1] == (-)) &&
                isa(root.args[3].args[2], Number)
                root.args[2] = root.args[1](root.args[2], root.args[3].args[2])
                root.args[1] = root.args[3].args[1]
                root.args[3] = root.args[3].args[3]
            end
            if (root.args[1] == (+) || root.args[1] == (-)) && # (x ± y) ± z = x ± (y ± z) = x ± w
                isa(root.args[3], Number) &&
                isa(root.args[2], Expr) && (root.args[2].args[1] == (+) || root.args[2].args[1] == (-)) &&
                isa(root.args[2].args[3], Number)
                root.args[3] = root.args[1](root.args[3], root.args[2].args[3])
                root.args[1] = root.args[2].args[1]
                root.args[2] = root.args[2].args[2]
            end
            if all(e->!(isa(e, Symbol) || isa(e, Expr)), root.args[2:3]) # arithmetic
                root = root.args[1](root.args[2:3]...)
            end
        end
    end
    return root
end

function infix(io::IO, root; digits=3)
    if isa(root, Number)
        print(io, round(root, digits=digits))
    elseif isa(root, Expr) && root.head == :call
        print(io, "(")
        infix(io, root.args[2])
        infix(io, root.args[1])
        infix(io, root.args[3])
        print(io, ")")
    else
        print(io, root)
    end
end

function latex(root; digits=3)
    if isa(root, Number)
        print(round(root, digits=digits))
    elseif isa(root, Expr) && root.head == :call
        if root.args[1] == (/)
            print("\\frac{")
            latex(root.args[2])
            print("}{")
            latex(root.args[3])
            print("}")
        else
            print("\\left(")
            latex(root.args[2])
            latex(root.args[1])
            latex(root.args[3])
            print("\\right)")
        end
    else
        print(root)
    end
end

# Custom optimization call
function optimize(f, method::TreeGP, options::Options = Options(;default_options(method)...))
    method.optimizer.mutation =  mutatetree(method)
    method.optimizer.crossover = crosstree
    optimize(f, NoConstraints(), nothing, method, options)
end


text:
mutable struct Env

    config :: Dict

    sub_yrng :: Union{UnitRange, Colon, Nothing}
    Nx :: Integer
    Ny :: Integer
    Nz :: Integer

    cdata_varnames :: AbstractArray{String}

    gf        :: Union{PolelikeCoordinate.GridFile, Nothing}
    gd        :: PolelikeCoordinate.Grid
    gd_slab   :: PolelikeCoordinate.Grid
    topo      :: Topography

    function Env(
        config;
        sub_yrng :: Union{UnitRange, Colon} = Colon(),
        verbose :: Bool = false,
    )

        writeLog("Validating config: MODEL_CORE")
        config = validateConfigEntries(config, getConfigDescriptor()["MODEL_CORE"]; verbose = verbose)

        # mask =>   lnd = 0, ocn = 1
        gf = PolelikeCoordinate.CurvilinearSphericalGridFile(
            config["domain_file"];
            R  = Re,
            Ω  = Ω,
        )

        if sub_yrng == Colon()
            sub_yrng = 1:gf.Ny
        end

        Nx = gf.Nx
        Ny = length(sub_yrng)
        Nz = length(config["z_w"]) - 1

        cdata_varnames = []

        if config["MLD_scheme"] == "datastream"
            push!(cdata_varnames, "HMXL")
        end

        if config["Qflx"] == "on"
            push!(cdata_varnames, "QFLXT")
            push!(cdata_varnames, "QFLXS")
        end

        if config["weak_restoring"] == "on" || config["Qflx_finding"] == "on"
            push!(cdata_varnames, "TEMP")
            push!(cdata_varnames, "SALT")
        end

        # I put gd here because master needs gd for output coordinate information
        gd      = PolelikeCoordinate.genGrid(gf, config["z_w"] ; sub_yrng=sub_yrng)
        gd_slab = PolelikeCoordinate.genGrid(gf, [0, -1.0]; sub_yrng=sub_yrng)

        if config["Returnflow_layers"] + config["Ekman_layers"] > gd.Nz
            throw(ErrorException("Sum of `Returnflow_layers` and `Ekman_layers` exceeds Nz = $(gd.Nz)"))
        end
        #
        topo = nothing

        if config["topo_file"] != ""
            Dataset(config["topo_file"], "r") do ds
                Nz_bot = ds["Nz_bot"][:, sub_yrng]
            end
        else
            Nz_bot = zeros(Int64, gd.Nx, gd.Ny) .+ gd.Nz
            Nz_bot[gf.mask[:, sub_yrng] .== 0.0 ] .= 0.0
        end
        topo = Topography(
            Nz_bot, gf.Nx, length(sub_yrng), config["z_w"];
            deep_depth = - config["z_w"][config["Returnflow_layers"] + config["Ekman_layers"] + 1]
        )

        if any(view(topo.mask_T, 1, :, :) .!= gf.mask[:, sub_yrng])
            throw(ErrorException("Topo surface mask is inconsistent with mask loaded from domain file. It can be that the `deep_depth` in Topography is set too low, or simply because topography file is not consistent with domain file."))
        end

        return new(

            config,

            sub_yrng,
            Nx,
            Ny,
            Nz,

            cdata_varnames,

            gf,
            gd,
            gd_slab,
            topo,
        )
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
export inv, transpose, translate, scale, has_scale,
rotate_x, rotate_y, rotate_z, rotate, look_at, swaps_handedness

Base.inv(t::Transform{T}) where T<:AbstractFloat = Transform{T}(t.m_inv, t.m)

Base.transpose(t::Transform{T}) where T<:AbstractFloat =
    Transform{T}(transpose(t.m), transpose(t.m_inv))

function translate(Δ::Vector{3,T}) where T<:AbstractFloat
    m = SMatrix{4,4,T}(
        1.0, 0.0, 0.0, 0.0,
        0.0, 1.0, 0.0, 0.0,
        0.0, 0.0, 1.0, 0.0,
        Δ.x, Δ.y, Δ.z, 1.0)
    m_inv = SMatrix{4,4,T}(
         1.0,  0.0, 0.0, 0.0,
         0.0,  1.0, 0.0, 0.0,
         0.0,  0.0, 1.0, 0.0,
        -Δ.x, -Δ.y, Δ.z, 1.0)
    Transform{T}(m, m_inv)
end

function scale(x::T, y::U, z::V) where {T<:AbstractFloat, U<:AbstractFloat, V<:AbstractFloat}
    m = SMatrix{4,4,promote_type(T,U,V)}(
        x,   0.0, 0.0, 0.0,
        0.0, y,   0.0, 0.0,
        0.0, 0.0, z,   0.0,
        0.0, 0.0, 0.0, 1.0)
    m_inv = SMatrix{4,4,promote_type(T,U,V)}(
         1.0/x, 0.0,   0.0,   0.0,
         0.0,   1.0/y, 0.0,   0.0,
         0.0,   0.0,   1.0/z, 0.0,
         0.0,   0.0,   0.0,   1.0)
    Transform{T}(m, m_inv)
end

function has_scale(t::Transform{T}) where T<:AbstractFloat
    @. axes = magnitude² |> t |> (Vector3{T}(1,0,0), Vector3{T}(0,1,0), Vector3{T}(0,0,1))
    all(axes .≈ 1.0)
end

function rotate_x(θ::T) where T<:AbstractFloat
    sinθ = (sin ∘ deg2rad)(θ)
    cosθ = (cos ∘ deg2rad)(θ)
    m = SMatrix{4,4,T}(
        1.0,  0.0,  0.0,  0.0,
        0.0,  cosθ, sinθ, 0.0,
        0.0, -sinθ, cosθ, 0.0,
        0.0,  0.0,  0.0,  1.0)
    Transform{T}(m, transpose(m))
end

function rotate_y(θ::T) where T<:AbstractFloat
    sinθ = (sin ∘ deg2rad)(θ)
    cosθ = (cos ∘ deg2rad)(θ)
    m = SMatrix{4,4,T}(
        cosθ, 0.0, -sinθ, 0.0,
        0.0,  1.0,  0.0,  0.0,
        sinθ, 0.0,  cosθ, 0.0,
        0.0,  0.0,  0.0,  1.0)
    Transform{T}(m, transpose(m))
end

function rotate_z(θ::T) where T<:AbstractFloat
    sinθ = (sin ∘ deg2rad)(θ)
    cosθ = (cos ∘ deg2rad)(θ)
    m = SMatrix{4,4,T}(
         cosθ, sinθ, 0.0, 0.0,
        -sinθ, cosθ, 0.0, 0.0,
         0.0,  0.0,  1.0, 0.0,
         0.0,  0.0,  0.0, 1.0)
    Transform{T}(m, transpose(m))
end

function rotate(θ::T, axis::Vector{3,U}) where {T<:AbstractFloat, U<:AbstractFloat}
    aₙ = normalize(axis)
    sinθ = (sin ∘ deg2rad)(θ)
    cosθ = (cos ∘ deg2rad)(θ)
    m = SMatrix{4,4,promote_type(T,U)}(
        a.x * a.x + (1 - a.x * a.x) * cosθ,
        a.x * a.y * (1 - cosθ) + a.z * sinθ,
        a.x * a.z * (1 - cosθ) - a.y * sinθ,
        0.0,
        a.x * a.y * (1 - cosθ) - a.z * sinθ,
        a.y * a.y + (1 - a.y * a.y) * cosθ,
        a.y * a.z * (1 - cosθ) + a.x * sinθ,
        0.0,
        a.x * a.z * (1 - cosθ) + a.y * sinθ,
        a.y * a.z * (1 - cosθ) - a.x * sinθ,
        a.z * a.z + (1 - a.z * a.z) * cosθ,
        0.0,
        0.0, 0.0, 0.0, 1.0)
    Transform{T}(m, transpose(m))
end

# I'm very sorry. I've learned about adding latex symbols and I can't stop. Please understand.
function look_at(📷::Point{3,T}, at::Point{3,U}, up::Vector{3,V}) where {T<:AbstractFloat, U<:AbstractFloat, V<:AbstractFloat}
    ↗ = normalize(at - eye)
    ← = normalize(normalize(up) × direction)
    ↑ = ↗ × ←
    m = SMatrix{4,4,promote_type(T,U,V)}(
        ←.x,  ←.y,   ←.z,  0.0,
        ↑.x,  ↑.y,   ↑.z,  0.0,
        ↗.x,  ↗.y,   ↗.z,  0.0,
        📷.x, 📷.y, 📷.z, 1.0)
    Transform(inv(m), m)
end

Base.:(*)(t::Transform{T}, u::Transform{U}) where T = Transform{promote_type(T,U)}(t.m * u.m, u.m_inv * t.m_inv)

swaps_handedness(t::Transform{T}) where T = det(t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,3)]) < zero(T)

function Transform(q::Quaternion{T}) where T<:AbstractFloat
    m = SMatrix{4,4,T}(
        1 - 2 * (q.v.y * q.v.y + q.v.z * q.v.z),
        2 * (q.v.x * q.v.y + q.w * q.v.z),
        2 * (q.v.x * q.v.z - q.w * q.v.y),
        0.0,
        2 * (q.v.x * q.v.y - q.w * q.v.z),
        1 - 2 * (q.v.x * q.v.x + q.v.z * q.v.z),
        2 * (q.v.y * q.v.z + q.w * q.v.x),
        0.0,
        2 * (q.v.x * q.v.z + q.w * q.v.y),
        2 * (q.v.y * q.v.z - q.w * q.v.x),
        1 - 2 * (q.v.x * q.v.x + q.v.y * q.v.y),
        0.0,
        0.0, 0.0, 0.0, 1.0)
    Transform{T}(m, transpose(m))
end


text:
#=

contains can_copy and tm_hash

We need to be able to tm_hash every variable for which can_copy returns true.
The code for can_copy is naturally based on Julia's Base.deepcopy.
We use the same code to compute hashes because
1. We only need to hash variables we can copy.
2. We want to be able to hash circular data structures. Base.hash can fail on those.
3. We want a hash with the property that hash(x) != hash(y) means that it is unlikely that x == y.
Base.hash does not have this property.

=#

"""
    can_copy(x)

Is supposed to return true if deepcopy works on x.
The code to test this is based on deepcopy, but it could get some strange case wrong.
"""
can_copy(x) = can_copy_and_hash(x)[1]

"""
    tm_hash(x)

Produces a hash of any variable for which `can_copy` returns true.
If variables `x` and `y` are different, they will probably have different hashes.
"""
tm_hash(x) = can_copy_and_hash(x)[2]

can_copy_and_hash(x) = can_copy_and_hash(x, IdDict(), zero(UInt64))
can_copy_and_hash(x::Union{Core.MethodInstance,Module,DataType,Method,GlobalRef,UnionAll,Task,Regex,Function,IO}, id, h) = (false, zero(UInt64))
can_copy_and_hash(x::Union{Symbol, Union, String}, id, h) = (true, hash(x,h))

function can_copy_and_hash(x, id, h)
    isbitstype(typeof(x)) && return (true, hash(x,h))
    haskey(id, x) && return (true, hash(objectid(x),h))
    id[x] = true
    nf = nfields(x)
    nf == 0 && return (true, hash(x, h))
    id[x] = true
    val = true
    h = hash(typeof(x), h)
    for i in 1:nf
        tf, h = can_copy_and_hash(getfield(x,i),id,h)
        if tf == false
            val = false
            break
        end
    end
    return val, h
end

function can_copy_and_hash(x::Union{Tuple,Core.SimpleVector,Array}, id, h)
    haskey(id, x) && return (true, hash(objectid(x), h))
    id[x] = true
    h = hash(typeof(x), h)
    if isa(x, Array)
        h = hash(size(x), h)
    end
    val = true
    for i in 1:length(x)
        tf, h = can_copy_and_hash(x[i],id,h)
        if tf == false
            val = false
            break
        end
    end
    return val, h
end

function can_copy_and_hash(x::Dict, id, h)
    haskey(id, x) && return (true, hash(objectid(x), h))
    id[x] = true

    h = hash(typeof(x), h)
    val = true
    for (k,v) in x
        tf, h = can_copy_and_hash((k,v),id,h)
        if tf == false
            val = false
            break
        end
    end
    return val, h

end

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


text: print("hai")

max_stars_repo_name: urbanDemon/test

max_issues_repo_name: urbanDemon/test

max_forks_repo_name: urbanDemon/test


text:
module Phi0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    [-4.42055412371134e-10*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 1.65770779639175e-10*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 1.51956548002578e-10*xi.^2.*(l1 - 0.0001).^2 + 1.20560567010309e-10*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 8.0373711340206e-11*xi.^2.*(l2 - 0.002).^2 + 1.13025531572165e-11*xi.^2.*(l3 - 0.0015).^2, 1.16009459695093e-10*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) + 4.35035473856599e-11*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) - 3.98782517701883e-11*xi.^2.*(l1 - 0.0001).^2 - 3.16389435532073e-11*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) - 2.10926290354715e-11*xi.^2.*(l2 - 0.002).^2 - 2.96615095811319e-12*xi.^2.*(l3 - 0.0015).^2, -5.58213595360825e-5*xi.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 3.58060567010309e-5*xi.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 2.27779478092784e-5*xi.*(l1 - 0.0001).^2 - 3.898125e-8*xi.*(l1 - 0.0001) - 1.46850837628866e-5*xi.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 1.77203608247422e-5*xi.*(l2 - 0.002).^2 + 2.835e-8*xi.*(l2 - 0.002) + 1.26579413659794e-5*xi.*(l3 - 0.0015).^2 + 1.063125e-8*xi.*(l3 - 0.0015)]
end
end


text:
#
# Example script to wrap libLLVM (LLVM C interface)
#

using Clang.cindex
using Clang.wrap_c

LLVM_VER = "3.3"
LLVM_BUILD_TYPE="Release+Asserts"

JULIA_ROOT=abspath(JULIA_HOME, "../../")
LLVM_PATH=joinpath(JULIA_ROOT, "deps/llvm-$LLVM_VER")
libLLVM_PATH = joinpath(LLVM_PATH, "include/llvm-c")
clanginc_path = joinpath(LLVM_PATH, "build/$LLVM_BUILD_TYPE/lib/clang/3.3/include")

libllvm_hdrs = map(x->joinpath(libLLVM_PATH, x),
                    split(readall(`find $libLLVM_PATH -name *.h`))
                  )

clang_includes = map(x::ASCIIString->joinpath(LLVM_PATH, x), [
    "build/$LLVM_BUILD_TYPE/lib/clang/3.3/include",
    "build/include/",
    "include"
    ])
clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

const wc = wrap_c.init(;
                        output_file = "libLLVM_h.jl",
                        common_file = "libLLVM_common.jl",
                        clang_includes = clang_includes,
                        clang_args = clang_extraargs,
                        header_library = x->:libllvm,
                        header_wrapped = (x,y)->contains(y, "LLVM") )

function wrap_libLLVM(wc::WrapContext, wrap_hdrs)
    wrap_c.wrap_c_headers(wc, convert(Vector{ASCIIString}, wrap_hdrs))
end

if !isinteractive()
    wrap_libLLVM(wc, libllvm_hdrs)
end


text:
module Liu

using ModelingToolkit
#include("../../channel_types.jl")

using ..NeuronBuilder
import ..channel_dynamics
include("channels.jl")

export channel_dynamics

end


text:
export LBFGS, update!, reset!

"""
`LBFGS(domainType::Type,dim_in::Tuple, M::Integer)`

`LBFGS(dim_in::Tuple, M::Integer)`

`LBFGS(x::AbstractArray, M::Integer)`

Construct a Limited-Memory BFGS `LinearOperator` with memory `M`. The memory of `LBFGS` can be updated using the function `update!`, where the current iteration variable and gradient (`x`, `grad`) and the previous ones (`x_prev` and `grad_prev`) are needed:

```
julia> L = LBFGS(Float64,(4,),5)
LBFGS  ℝ^4 -> ℝ^4

julia> update!(L,x,x_prev,grad,grad_prev); # update memory

julia> d = L*grad; # compute new direction

```

Use  `reset!(L)` to cancel the memory of the operator.

"""
mutable struct LBFGS{R, T <: AbstractArray, M, I <: Integer} <: LinearOperator
  currmem::I
  curridx::I
  s::T
  y::T
  s_M::Array{T, 1}
  y_M::Array{T, 1}
  ys_M::Array{R, 1}
  alphas::Array{R, 1}
  H::R
end

#default constructor

function LBFGS(x::T, M::I) where {T <: AbstractArray, I <: Integer}
  s_M = [zero(x) for i = 1:M]
  y_M = [zero(x) for i = 1:M]
  s = zero(x)
  y = zero(x)
  R = real(eltype(x))
  ys_M = zeros(R, M)
  alphas = zeros(R, M)
  LBFGS{R, T, M, I}(0, 0, s, y, s_M, y_M, ys_M, alphas, one(R))
end

"""
`update!(L::LBFGS, x, x_prex, grad, grad_prev)`

See the documentation for `LBFGS`.
"""
function update!(L::LBFGS{R, T, M, I}, x::T, x_prev::T, gradx::T, gradx_prev::T) where {R, T, M, I}
  L.s .= x .- x_prev
  L.y .= gradx .- gradx_prev
  ys = real(dot(L.s, L.y))
  if ys > 0
    L.curridx += 1
    if L.curridx > M L.curridx = 1 end
    L.currmem += 1
    if L.currmem > M L.currmem = M end
    L.ys_M[L.curridx] = ys
    L.s_M[L.curridx] .= L.s
    L.y_M[L.curridx] .= L.y
    yty = real(dot(L.y, L.y))
    L.H = ys/yty
  end
  return L
end

"""
`reset!(L::LBFGS)`

Cancels the memory of `L`.
"""
function reset!(L::LBFGS)
  L.currmem, L.curridx = 0, 0
  L.H = 1.0
end

# LBFGS operators are symmetric

mul!(x::T, L::AdjointOperator{LBFGS{R, T, M, I}}, y::T) where {R, T, M, I} = mul!(x, L.A, y)

# Two-loop recursion

function mul!(d::T, L::LBFGS{R, T, M, I}, gradx::T) where {R, T, M, I}
  d .= gradx
  idx = loop1!(d,L)
  d .*= L.H
  loop2!(d,idx,L)
  return d
end

function loop1!(d::T, L::LBFGS{R, T, M, I}) where {R, T, M, I}
  idx = L.curridx
  for i = 1:L.currmem
    L.alphas[idx] = real(dot(L.s_M[idx], d))/L.ys_M[idx]
    d .-= L.alphas[idx] .* L.y_M[idx]
    idx -= 1
    if idx == 0 idx = M end
  end
  return idx
end

function loop2!(d::T, idx, L::LBFGS{R, T, M, I}) where {R, T, M, I}
  for i = 1:L.currmem
    idx += 1
    if idx > M idx = 1 end
    beta = real(dot(L.y_M[idx], d))/L.ys_M[idx]
    d .+= (L.alphas[idx] - beta) .* L.s_M[idx]
  end
  return d
end

# Properties

domainType(L::LBFGS)   =   eltype(L.y_M[1])
domainType(L::LBFGS{R,T}) where { R, T <: ArrayPartition } = eltype.(L.y_M[1].x)
codomainType(L::LBFGS) = eltype(L.y_M[1])
codomainType(L::LBFGS{R,T}) where { R, T <: ArrayPartition } = eltype.(L.y_M[1].x)

size(A::LBFGS{R,T}) where { R, T <: ArrayPartition } = (size.(A.s.x), size.(A.s.x))
size(A::LBFGS) = (size(A.s), size(A.s))

fun_name(A::LBFGS) = "LBFGS"


text:
function poynting_vector(k0, beta, centers, points, Ez_inc, Hx_inc, Hy_inc)
    #assumes points lie outside all scattering disks
    Ns = size(centers,1)
	P = div(div(length(beta),Ns)-1,2)
	len = size(points,1)
	pyntg = Array{Float64}(undef, len, 2)
    Ez = copy(Ez_inc)
    for ip in 1:len
        ∂Ez∂x = 0.0im
        ∂Ez∂y = 0.0im
        for ic in 1:Ns
		    ind = (ic-1)*(2*P+1) .+ P + 1
			pt1 = points[ip,1] - centers[ic,1]
			pt2 = points[ip,2] - centers[ic,2]

			R = hypot(pt1, pt2)
			θ = atan(pt2, pt1)
			R == 0 && error("R == 0, center=$(centers[ic,:]),
                                point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip")

            Hₚ₋₂ = besselh(-1, 1, k0*R)
            Hₚ₋₁ = besselh(0, 1, k0*R)

            Ez[ip] += beta[ind]*Hₚ₋₁
            ∂H∂R = k0*Hₚ₋₂
            ∂Ez∂x += beta[ind]*∂H∂R*pt1/R
            ∂Ez∂y += beta[ind]*∂H∂R*pt2/R

            for p = 1:P
                Hₚ = (2*(p-1)/(k0*R))*Hₚ₋₁ - Hₚ₋₂
                ∂H∂R = k0*(Hₚ₋₁ - p/(k0*R)*Hₚ)

                Ez[ip] += beta[ind - p]*(-1)^p*Hₚ*exp(-1im*p*θ)
                ∂Ez∂x += beta[ind - p]*(-1)^p*exp(-1im*p*θ)*(∂H∂R*pt1/R + Hₚ*1im*(-p)*(-pt2)/R^2)
                ∂Ez∂y += beta[ind - p]*(-1)^p*exp(-1im*p*θ)*(∂H∂R*pt2/R + Hₚ*1im*(-p)*(pt1)/R^2)

                Ez[ip] += beta[ind + p]*Hₚ*exp(1im*p*θ)
                ∂Ez∂x += beta[ind + p]*exp(1im*p*θ)*(∂H∂R*pt1/R + Hₚ*1im*p*(-pt2)/R^2)
                ∂Ez∂y += beta[ind + p]*exp(1im*p*θ)*(∂H∂R*pt2/R + Hₚ*1im*p*(pt1)/R^2)
                Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ
			end
		end
        Hx = Hx_inc[ip] + (1/1im/k0/eta0)*∂Ez∂y
        Hy = Hy_inc[ip] + (-1/1im/k0/eta0)*∂Ez∂x
        pyntg[ip,1] = (-0.5)*real(Ez[ip]*conj(Hy))
        pyntg[ip,2] = 0.5*real(Ez[ip]*conj(Hx))
	end
    pyntg
end

function poynting_vector(k0, points, Ez_inc, Hx_inc, Hy_inc)
    #assumes points lie outside all scattering disks
    len = size(points,1)
	pyntg = Array{Float64}(undef, len, 2)
    for ip in 1:len
        pyntg[ip,1] = (-0.5)*real(Ez_inc[ip]*conj(Hy_inc[ip]))
        pyntg[ip,2] = 0.5*real(Ez_inc[ip]*conj(Hx_inc[ip]))
	end
    pyntg
end

function calc_power(k0::Array, kin::Array, P, sp, points, nhat, ui; opt = FMMoptions())
    #this assumes points are equidistant. correct result only after multiplying by arc length
    power = Array{Float64}(undef, length(k0))
    for i in eachindex(k0)
        power[i] = calc_power(k0[i], kin[i], P, sp, points, nhat, ui; opt=opt)
    end
    power
end

function calc_power(k0::Number, kin::Number, P, sp, points, nhat, ui; opt = FMMoptions())
    #this assumes points are equidistant. correct result only after multiplying by arc length
    Ez_inc = uinc(k0, points, ui)
    Hx_inc = hxinc(k0, points, ui)
    Hy_inc = hyinc(k0, points, ui)
    if sp == nothing
        pyntg = poynting_vector(k0, points, Ez_inc, Hx_inc, Hy_inc)
    else
        if opt.FMM
            result = solve_particle_scattering_FMM(k0, kin, P, sp, ui, opt;
                get_inner = false, verbose = false)
            if result[2].isconverged == false
                @warn("FMM process did not converge")
                return
            end
            beta = result[1]
        else
            beta = solve_particle_scattering(k0, kin, P, sp, ui;
                        get_inner = false, verbose = false)
        end
        pyntg = poynting_vector(k0, beta, sp.centers, points, Ez_inc, Hx_inc, Hy_inc)
    end
    power_quadrature(nhat, pyntg)
end

function power_quadrature(n, S)
    #trapezoidal rule
    len = size(S,1)
    (dot(n, S) - 0.5*(n[1,1]*S[1,1] + n[1,2]*S[1,2] +
        n[len,1]*S[len,1] + n[len,2]*S[len,2]))/(len - 1)
end

function rect_border(b, Nx, Ny)
    #clockwise from top
    points1 = [range(b[1], stop=b[2], length=Nx)     b[4]*ones(Nx)]
    points2 = [b[2]*ones(Ny)     range(b[4], stop=b[3], length=Ny)]
    points3 = [range(b[2], stop=b[1], length=Nx)     b[3]*ones(Nx)]
    points4 = [b[1]*ones(Ny)     range(b[3], stop=b[4], length=Ny)]

    n = [[zeros(Nx) ones(Nx)], [ones(Ny) zeros(Ny)], [zeros(Nx) -ones(Nx)], [-ones(Ny) zeros(Ny)]]
    [points1, points2, points3, points4], n
end


text:
# ---
# title: 1589. Maximum Sum Obtained of Any Permutation
# id: problem1589
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/maximum-sum-obtained-of-any-permutation/description/>
# hidden: true
# ---
#
# We have an array of integers, `nums`, and an array of `requests` where
# `requests[i] = [starti, endi]`. The `ith` request asks for the sum of
# `nums[starti] + nums[starti + 1] + ... + nums[endi - 1] + nums[endi]`. Both
# `starti` and `endi` are _0-indexed_.
#
# Return _the maximum total sum of all requests **among all permutations** of_
# `nums`.
#
# Since the answer may be too large, return it **modulo** `109 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4,5], requests = [[1,3],[0,1]]
#     Output: 19
#     Explanation: One permutation of nums is [2,1,3,4,5] with the following result:
#     requests[0] -> nums[1] + nums[2] + nums[3] = 1 + 3 + 4 = 8
#     requests[1] -> nums[0] + nums[1] = 2 + 1 = 3
#     Total sum: 8 + 3 = 11.
#     A permutation with a higher total sum is [3,5,4,2,1] with the following result:
#     requests[0] -> nums[1] + nums[2] + nums[3] = 5 + 4 + 2 = 11
#     requests[1] -> nums[0] + nums[1] = 3 + 5  = 8
#     Total sum: 11 + 8 = 19, which is the best that you can do.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4,5,6], requests = [[0,1]]
#     Output: 11
#     Explanation: A permutation with the max total sum is [6,5,4,3,2,1] with request sums [11].
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3,4,5,10], requests = [[0,2],[1,3],[1,1]]
#     Output: 47
#     Explanation: A permutation with the max total sum is [4,10,5,3,2,1] with request sums [19,18,10].
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `1 <= n <= 105`
#   * `0 <= nums[i] <= 105`
#   * `1 <= requests.length <= 105`
#   * `requests[i].length == 2`
#   * `0 <= starti <= endi < n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#
#   Purpose: Demonstrates how to solve the problem
#
#     maximize x^0.2*y^0.8 + z^0.4 - x
#           st x + y + 0.5z = 2
#              x,y,z >= 0
#

using Mosek
printstream(msg::AbstractString) = print(msg)

csub = [ 4, 5, 1 ]
cval = [ 1.0, 1.0, -1.0]
asub = [ 1, 2, 3]
aval = [ 1.0, 1.0, 0.5]
numvar = 6
numcon = 1

# Create a task
maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))

    appendcons(task,numcon)
    appendvars(task,numvar)

    # Set up the linear part of the problem
    putclist(task,csub, cval)
    putarow(task,1, asub, aval)
    putvarboundslice(task,1, 6,
                     [ MSK_BK_FR, MSK_BK_FR, MSK_BK_FR, MSK_BK_FR, MSK_BK_FR, MSK_BK_FX ],
                     [ -Inf,      -Inf,      -Inf,      -Inf,      -Inf,      1.0],
                     [  Inf,       Inf,       Inf,       Inf,       Inf,      1.0])
    putconbound(task,1, MSK_BK_FX, 2.0, 2.0)

    # Input the cones
    appendcone(task,MSK_CT_PPOW, 0.2, [1, 2, 4])
    appendcone(task,MSK_CT_PPOW, 0.4, [3, 6, 5])

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Optimize the task
    optimize(task,"mosek://solve.mosek.com:30080")
    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta == MSK_SOL_STA_OPTIMAL
        # Output a solution
        x = getxx(task,MSK_SOL_ITR)
        println("Optimal solution: $x")
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
using Memoize
@memoize josephus(n::Integer, k::Integer, m::Integer=1) = n == m ? collect(0:m .- 1) : mod.(josephus(n - 1, k, m) + k, n)

@show josephus(41, 3)
@show josephus(41, 3, 5)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 319ae2a8-b501-4078-b183-4832364c02fa
begin
	using Pkg
	Pkg.add("LinearAlgebra")
	Pkg.add("ProximalOperators")
	Pkg.add("FileIO")
	Pkg.add("ImageIO")
	Pkg.add("CUDA")
	Pkg.add("Flux")
	Pkg.add("Zygote")
	Pkg.add(url="https://github.com/JuliaTomo/TomoForward.jl")
	Pkg.add(url="https://github.com/JuliaTomo/DeepReconstruction.jl")
end

# ╔═╡ b94b20bc-9250-4889-83f1-eb5daab12e9b
begin
	using LinearAlgebra
	using ProximalOperators
	using Base.Iterators
	using Printf
	using CUDA, Flux, Zygote
	using TomoForward, DeepReconstruction
	using FileIO, ImageIO
end

# ╔═╡ 03d3ac83-1be7-4ab4-892f-091cbd758c95
md"""
This is a Pluto notebook for trying things out
"""

# ╔═╡ 8520fa00-9b21-11ec-08bd-a1067d39ea44
# import packages
#begin
#	using LinearAlgebra
#	using ProximalOperators
#	using Base.Iterators
#	using Printf
#end

# ╔═╡ d0ef7765-c56c-41a4-9087-88cde949fc5d
# set variables for testing
begin
	T = Float32
	A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
	b = T[1.0, 2.0, 3.0, 4.0]
	m, n = size(A)
	R = real(T)
	lam = R(0.1) * norm(A' * b, Inf)
	x0 = zeros(R, n)
	y0 = similar(x0)
	f = LeastSquares(A, b)
	g = NormL1(lam)
	gamma=R(10) / opnorm(A)^2
	proxf!(xhat,uhat) = prox!(xhat,f,uhat,gamma)
	denoiser!(yhat,rhat) = prox!(yhat,g,rhat,gamma)
end

# ╔═╡ 6e85550d-349a-48eb-802c-b3809db16dd2
begin
	# define the pnp-DRS interation
	Base.@kwdef struct PnpDrsIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},F<:Function,G<:Function}
	    proxf!::F
	    denoiser!::G
	    uhat0::Tx #uhat0 takes the place of x0
	    gamma::R #don't think gamma is needed
	end

	Base.IteratorSize(::Type{<:PnpDrsIteration}) = Base.IsInfinite()

	# define the pnp-DRS state
	Base.@kwdef struct PnpDrsState{Tx}
	    xhat::Tx
	    yhat::Tx = similar(xhat)
	    rhat::Tx = similar(xhat)
	    uhat::Tx = similar(xhat)
	end

	function Base.iterate(iter::PnpDrsIteration, state::PnpDrsState = PnpDrsState(xhat=copy(iter.uhat0)))
	    iter.proxf!(state.xhat,state.uhat) # xhat^{k+1} = A(uhat^{k+1})
	    state.rhat .= 2 .* state.xhat .- state.uhat
	    iter.denoiser!(state.yhat,state.rhat) # yhat^{k+1} = B(2xhat^{k+1} - uhat^k)
	    state.uhat .-= (state.xhat - state.yhat) #uhat^{k+1} = uhat^k - (xhat^{k+1} - yhat^{k+1})
	    return state, state
	end
end

# ╔═╡ 32c39089-40cf-49ce-bb20-5efbf759a0b5
#example: create an instance
pnp_dr_iter = PnpDrsIteration(proxf! = proxf!, denoiser! = denoiser!,uhat0=x0,gamma=gamma)

# ╔═╡ 22027231-0c43-4bd1-a21a-4b9f4a333568
begin
	img = Float32.(load("2foam30.png"))
	img_gt = copy(img)
end

# ╔═╡ ab0460df-6e85-4cc4-a4b1-3a37dbc05f59
begin
	H, W = size(img)
	detcount = H
	nangles = 30
	proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])
end

# ╔═╡ fd079dc5-403e-4a97-af12-60c489ad03f2
begin
	B = fp_op_parallel2d_line(proj_geom, size(img, 1), size(img, 2))
	p_data = Float32.(B * vec(img))
	# p_data = vec(reshape(p_data, nangles, :)') # detector count should be the first axis
	lr = 0.0001f0
end

# ╔═╡ 7607323c-a19e-48ab-9ec9-8a974c9a22e9
begin
	net = UNet(3; use_two_skip=true) # 128, 5
	opt = ADAM(lr) # Gadelha
	dresult = "../result/"
	mkpath(dresult)
	u_out, losses, u_best, errs = recon2d_dip(net, opt, p_data, Float32.(B), H, W; img_gt=img_gt)
end

# ╔═╡ c31ecec5-345a-4f3d-8653-410c908ca4bd
typeof(Float32.(B))

# ╔═╡ Cell order:
# ╟─03d3ac83-1be7-4ab4-892f-091cbd758c95
# ╠═319ae2a8-b501-4078-b183-4832364c02fa
# ╠═b94b20bc-9250-4889-83f1-eb5daab12e9b
# ╠═8520fa00-9b21-11ec-08bd-a1067d39ea44
# ╠═d0ef7765-c56c-41a4-9087-88cde949fc5d
# ╠═6e85550d-349a-48eb-802c-b3809db16dd2
# ╠═32c39089-40cf-49ce-bb20-5efbf759a0b5
# ╠═22027231-0c43-4bd1-a21a-4b9f4a333568
# ╠═ab0460df-6e85-4cc4-a4b1-3a37dbc05f59
# ╠═fd079dc5-403e-4a97-af12-60c489ad03f2
# ╠═7607323c-a19e-48ab-9ec9-8a974c9a22e9
# ╠═c31ecec5-345a-4f3d-8653-410c908ca4bd


text:
module PlotUtils

path_to_imports = joinpath(@__DIR__, "imports.jl")
include(path_to_imports)

# Load python defs
path_to_plot_defs = @__DIR__
pushfirst!(PyCall.PyVector(PyCall.pyimport("sys")."path"), "$path_to_plot_defs")
plot_yz = PyCall.pyimport("plot_yz")
blue2red = PyCall.pyimport("blue2red")
pyrange(n) = collect(range(0, stop=n-1))

# General plot settings
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 15
rcParams["xtick.labelsize"] = 15
rcParams["ytick.labelsize"] = 15
rcParams["figure.figsize"] = (6, 6)

extract(s, out) = [o[s] for o in out]

function boxplot(x; showmeans=true, whis=[2.5, 97.5], showfliers=false, kw...)
  plt.boxplot(x, showmeans=showmeans, whis=whis, showfliers=showfliers; kw...)
end

function add_gridlines_Z(Z)
  J, K = size(Z)
  for j in pyrange(J)
    plt.axhline(y=j+.5, color="grey", linewidth=.5)
  end

  for k in pyrange(K)
    plt.axvline(x=k+.5, color="grey", linewidth=.5)
  end
end

axhlines(x; kw...) = for xi in x plt.axhline(xi; kw...) end
cm_greys = plt.cm.get_cmap("Greys", 5)

# TODO
# Steal things from: Cytof5/sims/repfam_fs/test

function plot_Z(Z; colorbar=true)
  J, K = size(Z)
  p = plt.imshow(Z, aspect="auto", vmin=0, vmax=1, cmap=cm_greys)
  add_gridlines_Z(Z)
  plt.yticks(pyrange(J), pyrange(J) .+ 1, fontsize=rcParams["font.size"])
  plt.xticks(pyrange(K), pyrange(K) .+ 1, fontsize=rcParams["font.size"],
             rotation=90)
  if colorbar
    plt.colorbar()
  end
  return p
end

# plot y/z
function make_yz(y, Zs, Ws, lams, imgdir; vlim,
                 w_thresh=.01, lw=3,
                 Z_true=nothing,
                 markernames=[],
                 fs_y=rcParams["font.size"],
                 fs_z=rcParams["font.size"],
                 fs_ycbar=rcParams["font.size"],
                 fs_zcbar=rcParams["font.size"])
  # Make img dir if needed
  mkpath(imgdir)

  I = length(y)
  for i in 1:I
    idx_best = estimate_ZWi_index(Zs, Ws, i)

    Zi = Int.(Zs[idx_best])
    Wi = Float64.(Ws[idx_best][i, :])
    lami = Int64.(lams[idx_best][i])

    yi = Float64.(y[i])

    # plot Yi, lami
    plt.figure(figsize=(6, 6))
    plot_yz.plot_y(yi, Wi, lami, vlim=vlim, cm=blue2red.cm(9), lw=lw,
                   fs_xlab=fs_y, fs_ylab=fs_y, fs_lab=fs_y, fs_cbar=fs_ycbar,
                   markernames=markernames)
    plt.savefig("$(imgdir)/y$(i).pdf", bbox_inches="tight")
    plt.close()

    # plot Zi, Wi
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z(Zi, Wi, lami, w_thresh=w_thresh, add_colorbar=false,
                   fs_lab=fs_z, fs_celltypes=fs_z, fs_markers=fs_z,
                   fs_cbar=fs_zcbar, markernames=markernames)
    plt.savefig("$(imgdir)/Z$(i).pdf", bbox_inches="tight")
    plt.close()
  end

  if Z_true != nothing
    # plot Z true
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z_only(Z_true, fs=fs_z,
                        xlab="cell subpopulations", ylab="markers",
                        rotate_xticks=false)
    plt.savefig("$(imgdir)/Z_true.pdf", bbox_inches="tight")
    plt.close()

    # plot ZT true
    plt.figure(figsize=(6, 6))
    plot_yz.plot_Z_only(Z_true', fs=fs_z,
                        xlab="markers", ylab="cell subpopulations")
    plt.savefig("$(imgdir)/ZT_true.pdf", bbox_inches="tight")
    plt.close()
  end
end

grepKmcmc(s) = match(r"(?<=Kmcmc_)\d+", s)

function make_metrics(different_K_runs_dir, outputfname; thresh,
                      grepKmcmcfn=grepKmcmc)
  metrics = Dict{Int, Any}()

  for (root, dirs, files) in walkdir(different_K_runs_dir)
    for file in files
      if file == outputfname
        # Get path to output
        path_to_output = joinpath(root, file)
        println(path_to_output)

        # Parse Kmcmc
        Kmcmc = parse(Int, grepKmcmcfn(splitdir(root)[2]).match)

        # Load output
        output = BSON.load(path_to_output)

        # Append metrics to metrics dictionary
        metrics[Kmcmc] = output[:metrics]

        # Calibration metric
        cmetric = let
          Wmean = mean(extract(:W, output[:out][1]))
          num_wik_lt_thresh = sum(Wmean .< thresh)
        end

        # Append calibration metric to metrics dictionary
        metrics[Kmcmc][:cmetric] = cmetric

      end
    end
  end

  # All Kmcmcs
  Kmcmcs = sort(collect(keys(metrics)))

  # Create dir for metrics output
  metrics_dir = joinpath(different_K_runs_dir, "metrics")
  mkpath(metrics_dir)

  # Plot Calibration metric
  plt.plot([metrics[K][:cmetric] for K in Kmcmcs],
           [metrics[K][:LPML] for K in Kmcmcs], marker="o")
  for K in Kmcmcs
    plt.text(metrics[K][:cmetric], metrics[K][:LPML], K)
  end
  plt.xlabel(L"number of $W_{i,k}$ < 10%")
  plt.ylabel("LPML")
  plt.savefig(joinpath(metrics_dir, "calibration_metric.pdf"),
              bbox_inches="tight")
  plt.close()

  # Plot LPML, DIC
  for metric in [:LPML, :DIC]
    plt.plot(Kmcmcs, map(K -> metrics[K][metric], Kmcmcs), marker="o")
    plt.xlabel("K")
    plt.ylabel(String(metric))
    plt.savefig(joinpath(metrics_dir, "$(metric).pdf"), bbox_inches="tight")
    plt.close()
  end

  return metrics
end

end  # PlotUtils

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
## Gather data back from RDA files
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_1/kenya_inc_crit_no_vac_w_vac.rda")
"""
kenya_inc_crit_sc_1 = @rget kenya_inc_crit
kenya_inc_crit_sc_1 = sum(kenya_inc_crit_sc_1, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_2/kenya_inc_crit_30_perc_w_vac.rda")
"""
kenya_inc_crit_sc_2 = @rget kenya_inc_crit
kenya_inc_crit_sc_2 = sum(kenya_inc_crit_sc_2, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_3/kenya_inc_crit_50_perc_w_vac.rda")
"""
kenya_inc_crit_sc_3 = @rget kenya_inc_crit
kenya_inc_crit_sc_3 = sum(kenya_inc_crit_sc_3, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_4/kenya_inc_crit_70_perc_w_vac.rda")
"""
kenya_inc_crit_sc_4 = @rget kenya_inc_crit
kenya_inc_crit_sc_4 = sum(kenya_inc_crit_sc_4, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_5/kenya_inc_crit_30_perc_rapid_w_vac.rda")
"""
kenya_inc_crit_sc_5 = @rget kenya_inc_crit
kenya_inc_crit_sc_5 = sum(kenya_inc_crit_sc_5, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_6/kenya_inc_crit_50_perc_rapid_w_vac.rda")
"""
kenya_inc_crit_sc_6 = @rget kenya_inc_crit
kenya_inc_crit_sc_6 = sum(kenya_inc_crit_sc_6, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_7/kenya_inc_crit_70_perc_rapid_w_vac.rda")
"""
kenya_inc_crit_sc_7 = @rget kenya_inc_crit
kenya_inc_crit_sc_7 = sum(kenya_inc_crit_sc_7, dims = 2)[:, 1, :]
kenya_inc_crit = 0
R"""
kenya_inc_crit = 0 #Free memory
"""


text:
#!/bin/julia

using JSON3
using Base.Iterators: flatten

@enum GENDER MASCULINE FEMININE NEUTER

french = JSON3.read(read("$(@__DIR__)/../data/kaikki.org-dictionary-French.json"); jsonlines=true)

nouns =  filter(w -> w.pos == "noun" && !occursin(' ', w.word), french)

function sense2mf(tags)
    if "masculine" in tags
        MASCULINE
    elseif "feminine" in tags
        FEMININE
    else
        NEUTER
    end
end

sensegender = [[(word = n.word, gender = sense2mf(s.tags)) for s in n.senses if :tags in keys(s)] for n in nouns] |> flatten |> collect

# Last two characters of each word
Dict(s.word[prevind(s.word, ncodeunits(s.word), 1):end] => s.gender for s in sensegender if length(s.word) >= 3)

# Probably want to remove all the accents first:
# (In JS, do it this way)
#
#
# https://stackoverflow.com/questions/990904/remove-accents-diacritics-in-a-string-in-javascript#answers-header
# const str = "Crème Brulée"
# str.normalize("NFD").replace(/[\u0300-\u036f]/g, "")
# > "Creme Brulee"

# Fun things to do:
#
# Identify common endings, e.g. up to 5 characters long
#
# Plot gender versus ending - look for patterns


text:
using Base.Test
using Gunrock

a = include("init.jl")

#BFS
@test bfs(a,1) == [1, 2, 2, 2, 3, 3, 3]

#Node betweenness centrality
@test bc(a,1) == [0.0, 0.25, 0.5, 0.75, 0.0, 0.0, 0.0]

#Connected components
num_cc, ids = cc(a)
@test num_cc == 1
@test ids == [1, 1, 1, 1, 1, 1, 1]

#Single source shortest paths
@test sssp(a,1) == [0, 1, 1, 1, 2, 2, 2]

#Pagerank
nodes, ranks = pagerank(a)
@test nodes == [3, 4, 5, 6, 1, 2, 7]
@test_approx_eq ranks [1.2522, 1.0280, 1.0280, 1.0201, 0.8048, 0.8048, 0.7991]


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

# Variogram functions using local anisotropies

function qmat(q,m)
  # Anisotropy matrix for Mahalanobis distance
  N = length(m)
  P = quat_to_dcm(q)[SOneTo(N),SOneTo(N)]
  Λ = Diagonal(SVector{N}(one(eltype(m))./m.^2))
  Q = P'*Λ*P
  # need to explicitly set Symmetric(Q), but then det(Q) will fail with StaticArrays
  # while not fixed, setting Symmetric only when calling Mahalanobis
  # https://github.com/JuliaArrays/StaticArrays.jl/issues/955
end

function mwvario(estimator, localpar)
  # get local Mahalanobis matrix
  Q = qmat(localpar[1],localpar[2])
  # get reference pars. and apply local anisotropy to given structures
  p = structures(estimator.γ)
  γs = map(p[3]) do γ
    Qs = Q ./ radii(γ.ball)' .^ 2
    γ = @set γ.ball = MetricBall(1.0, Mahalanobis(Symmetric(Qs)))
  end
  γl = NuggetEffect(p[1]) + sum(c*γ for (c, γ) in zip(p[2], γs))

  # return local estimator
  if typeof(estimator) <: SimpleKriging
    return SimpleKriging(γl, estimator.mean)
  elseif typeof(estimator) <: OrdinaryKriging
    return OrdinaryKriging(γl)
  end

end

function kcfill!(Γ, γ::Variogram, X, localaniso)
  n = length(X)
  @inbounds for j in 1:n
    xj = X[j]
    Qj = localaniso[j]
    for i=j+1:n
      xi = X[i]
      Qi = localaniso[i]
      Γ[i,j] = kccov(γ, xi, xj, Qi, Qj)
    end
    Γ[j,j] = sill(γ)
    for i=1:j-1
      Γ[i,j] = Γ[j,i]
    end
  end
end

function kccov(γ::Variogram, xi, xj, Qi::AbstractMatrix, Qj::AbstractMatrix)
  Qij = (Qi+Qj)/2

  ## modify variogram with average anisotropy matrix
  p = structures(γ)
  γs = map(p[3]) do γx
    Qs = Qij ./ radii(γx.ball)' .^ 2
    γx = @set γx.ball = MetricBall(1.0, Mahalanobis(Symmetric(Qs)))
  end
  γl = NuggetEffect(p[1]) + sum(c*γx for (c, γx) in zip(p[2], γs))

  Cij = (sill(γl)-γl(xi,xj))
  (det(Qi)^0.25)*(det(Qj)^0.25)*(det(Qij)^-0.5)*Cij
end


text:
# Methods to provide grids and distance-matrices for
# different topologies of Kohonen maps (SOMs).

"""
    gridRectangular(xdim, ydim)

Create coordinates of all neurons on a rectangular SOM.

The return-value is an array of size (Number-of-neurons, 2) with
x- and y- coordinates of the neurons in the first and second
column respectively.
The distance between neighbours is 1.0.
The point of origin is bottom-left.
The first neuron sits at (0,0).
"""
function gridRectangular(xdim, ydim)

    grid = zeros(Float64, (xdim*ydim, 2))
    for ix in 1:xdim
        for iy in 1:ydim

            grid[ix+(iy-1)*xdim, 1] = ix-1
            grid[ix+(iy-1)*xdim, 2] = iy-1
        end
    end
    return grid
end

"""
    gridHexagonal(xdim, ydim)

Create coordinates of all neurons on a hexagonal SOM.

The return-value is an array of size (Number-of-neurons, 2) with
x- and y- coordinates of the neurons in the first and second
column respectively.
The distance between neighbours is 1.0.
The point of origin is bottom-left.
The first neuron sits at (0,0).
"""
function gridHexagonal(xdim, ydim)

    grid = zeros(Float64, (xdim*ydim, 2))
    rowHeight = √(1-0.25)

    for ix in 1:xdim
        for iy in 1:ydim

            if iy % 2 != 0    # 1st, 3rd, ... line
                grid[ix+(iy-1)*xdim, 1] = ix - 1
            else                    # 2nd, 4th, ... line
                grid[ix+(iy-1)*xdim, 1] = ix - 0.5
            end

            grid[ix+(iy-1)*xdim, 2] = (iy-1) * rowHeight
        end
    end
    return grid
end

"""
    gridSpherical(numNeurons)

Create coordinates of dim neurons on a spherical SOM.

The neurons are almost equally distributed on a sphere by
applying Fibonacci's spirals.
The radius of the sphere is adapted to get a distance of neighbours
of about 1.0.
The return-value is an array of size (Number-of-neurons, 3) with
x-, y- and z- coordinates of the neurons in the first, second and third
column respectively.
The point of origin is bottom-left.
The center of the sphere sits at (0,0).
"""
function gridSpherical(numNeurons)

    r = √(numNeurons/(4 * π))
    grid = zeros(Float64, (numNeurons, 3))

    # iterate Fibonacci spiral:
    for i in 1:(numNeurons)
        py = (i-1) * 2/numNeurons - 1 + (1/numNeurons)
        pr = √(1-py^2)
        ϕ  = (i-1) * π * (3-√5)
        px = cos(ϕ) * pr
        pz = sin(ϕ) * pr

        #println("raw: $px $py $pz")

        # shift point from origin to atom location:
        px = r * px
        py = r * py
        pz = r * pz

        #println("scaled: $px $py $pz")
        grid[i,:] = [px, py, pz]
    end
    return grid
end

# verify sphere:
# x = gridSpherical(50)
# R"""
# library(rgl)
# plot3d($x)
# #readline(prompt="Press [enter] to continue")
# """

"""
    distMatrix(grid::Array, toroidal::Bool)

Return the distance matrix for a non-toroidal or toroidal SOM.

# Arguments
- `grid`: coordinates of all neurons as generated by
          one of the `grid-`functions with x-coordinates in 1st column
          and y-coordinates in 2nd column.
- `toroidal`: true for a toroidal SOM.
"""
function distMatrix(grid::Array, toroidal::Bool)

    X = 1
    Y = 2
    xdim = maximum(grid[:,X]) - minimum(grid[:,X]) + 1.0
    ydim = maximum(grid[:,Y]) - minimum(grid[:,Y]) + 1.0

    numNeurons = nrow(grid)

    dm = zeros(Float64, (numNeurons,numNeurons))
    for i in 1:numNeurons
        for j in 1:numNeurons
            Δx = abs(grid[i,X] - grid[j,X])
            Δy = abs(grid[i,Y] - grid[j,Y])

            if toroidal
                Δx = min(Δx, xdim-Δx)
                Δy = min(Δy, ydim-Δy)
            end

            dm[i,j] = √(Δx^2 + Δy^2)
        end
    end
    # show(STDOUT, "text/plain",  grid)
    # show(STDOUT, "text/plain",  dm)
    return dm
end

"""
    distMatrixSphere(grid::Array)

Return the distance matrix for a spherical SOM.

# Arguments
- `grid`: coordinates of all neurons as gereated by
          one of the `grid-`functions with x-coordinates in 1st,
          y-coordinates in 2nd and z-coordinates in 3rd column.
          Point of origin must be in the center of the spere.
"""
function distMatrixSphere(grid::Array)

    r = (maximum(grid) - minimum(grid)) / 2

    numNeurons = nrow(grid)

    dm = zeros(Float64, (numNeurons,numNeurons))
    for i in 1:numNeurons
        for j in 1:numNeurons
            len_i = norm(grid[i,:])
            len_j = norm(grid[j,:])
            bog = grid[i,:] ⋅ grid[j,:] / (len_i * len_j)
            if bog > 1
                bog = 1.0
            end
            α = acos(bog)
            dm[i,j] = r * α
        end
    end
    return dm
end