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

# These tests cover the higher order functions specialized for sparse arrays defined in
# base/sparse/higherorderfns.jl, particularly map[!]/broadcast[!] for SparseVectors and
# SparseMatrixCSCs at present.

@testset "map[!] implementation specialized for a single (input) sparse vector/matrix" begin
    N, M = 10, 12
    for shapeA in ((N,), (N, M))
        A = sprand(shapeA..., 0.4); fA = Array(A)
        # --> test map entry point
        @test map(sin, A) == sparse(map(sin, fA))
        @test map(cos, A) == sparse(map(cos, fA))
        # --> test map! entry point
        fX = copy(fA); X = sparse(fX)
        map!(sin, X, A); X = sparse(fX) # warmup for @allocated
        @test (@allocated map!(sin, X, A)) == 0
        @test map!(sin, X, A) == sparse(map!(sin, fX, fA))
        @test map!(cos, X, A) == sparse(map!(cos, fX, fA))
        @test_throws DimensionMismatch map!(sin, X, spzeros((shapeA .- 1)...))
    end
end

@testset "map[!] implementation specialized for a pair of (input) sparse vectors/matrices" begin
    N, M = 10, 12
    f(x, y) = x + y + 1
    for shapeA in ((N,), (N, M))
        A, Bo = sprand(shapeA..., 0.3), sprand(shapeA..., 0.3)
        B = ndims(Bo) == 1 ? SparseVector{Float32, Int32}(Bo) : SparseMatrixCSC{Float32,Int32}(Bo)
        # use different types to check internal type stability via allocation tests below
        fA, fB = map(Array, (A, B))
        # --> test map entry point
        @test map(+, A, B) == sparse(map(+, fA, fB))
        @test map(*, A, B) == sparse(map(*, fA, fB))
        @test map(f, A, B) == sparse(map(f, fA, fB))
        @test_throws DimensionMismatch map(+, A, spzeros((shapeA .- 1)...))
        # --> test map! entry point
        fX = map(+, fA, fB); X = sparse(fX)
        map!(+, X, A, B); X = sparse(fX) # warmup for @allocated
        @test (@allocated map!(+, X, A, B)) == 0
        @test map!(+, X, A, B) == sparse(map!(+, fX, fA, fB))
        fX = map(*, fA, fB); X = sparse(fX)
        map!(*, X, A, B); X = sparse(fX) # warmup for @allocated
        @test (@allocated map!(*, X, A, B)) == 0
        @test map!(*, X, A, B) == sparse(map!(*, fX, fA, fB))
        @test map!(f, X, A, B) == sparse(map!(f, fX, fA, fB))
        @test_throws DimensionMismatch map!(f, X, A, spzeros((shapeA .- 1)...))
    end
end

@testset "map[!] implementation capable of handling >2 (input) sparse vectors/matrices" begin
    N, M = 10, 12
    f(x, y, z) = x + y + z + 1
    for shapeA in ((N,), (N, M))
        A, B, Co = sprand(shapeA..., 0.2), sprand(shapeA..., 0.2), sprand(shapeA..., 0.2)
        C = ndims(Co) == 1 ? SparseVector{Float32,Int32}(Co) : SparseMatrixCSC{Float32,Int32}(Co)
        # use different types to check internal type stability via allocation tests below
        fA, fB, fC = map(Array, (A, B, C))
        # --> test map entry point
        @test map(+, A, B, C) == sparse(map(+, fA, fB, fC))
        @test map(*, A, B, C) == sparse(map(*, fA, fB, fC))
        @test map(f, A, B, C) == sparse(map(f, fA, fB, fC))
        @test_throws DimensionMismatch map(+, A, B, spzeros(N, M - 1))
        # --> test map! entry point
        fX = map(+, fA, fB, fC); X = sparse(fX)
        map!(+, X, A, B, C); X = sparse(fX) # warmup for @allocated
        @test (@allocated map!(+, X, A, B, C)) == 0
        @test map!(+, X, A, B, C) == sparse(map!(+, fX, fA, fB, fC))
        fX = map(*, fA, fB, fC); X = sparse(fX)
        map!(*, X, A, B, C); X = sparse(fX) # warmup for @allocated
        @test (@allocated map!(*, X, A, B, C)) == 0
        @test map!(*, X, A, B, C) == sparse(map!(*, fX, fA, fB, fC))
        @test map!(f, X, A, B, C) == sparse(map!(f, fX, fA, fB, fC))
        @test_throws DimensionMismatch map!(f, X, A, B, spzeros((shapeA .- 1)...))
    end
end

@testset "broadcast! implementation specialized for solely an output sparse vector/matrix (no inputs)" begin
    N, M, p = 10, 12, 0.4
    V, C = sprand(N, p), sprand(N, M, p)
    fV, fC = Array(V), Array(C)
    @test broadcast!(() -> 0, V) == sparse(broadcast!(() -> 0, fV))
    @test broadcast!(() -> 0, C) == sparse(broadcast!(() -> 0, fC))
    @test let z = 0, fz = 0; broadcast!(() -> z += 1, V) == broadcast!(() -> fz += 1, fV); end
    @test let z = 0, fz = 0; broadcast!(() -> z += 1, C) == broadcast!(() -> fz += 1, fC); end
end

@testset "broadcast implementation specialized for a single (input) sparse vector/matrix" begin
    # broadcast for a single (input) sparse vector/matrix falls back to map, tested
    # extensively above. here we simply lightly exercise the relevant broadcast entry
    # point.
    N, M, p = 10, 12, 0.4
    a, A = sprand(N, p), sprand(N, M, p)
    fa, fA = Array(a), Array(A)
    @test broadcast(sin, a) == sparse(broadcast(sin, fa))
    @test broadcast(sin, A) == sparse(broadcast(sin, fA))
end

@testset "broadcast! implementation specialized for a single (input) sparse vector/matrix" begin
    N, M, p = 10, 12, 0.3
    f(x, y) = x + y + 1
    mats = (sprand(N, M, p), sprand(N, 1, p), sprand(1, M, p), sprand(1, 1, 1.0), spzeros(1, 1))
    vecs = (sprand(N, p), sprand(1, 1.0), spzeros(1))
    # --> test with matrix destination (Z/fZ)
    fZ = Array(first(mats))
    for Xo in (mats..., vecs...)
        X = ndims(Xo) == 1 ? SparseVector{Float32,Int32}(Xo) : SparseMatrixCSC{Float32,Int32}(Xo)
        shapeX, fX = size(X), Array(X)
        # --> test broadcast! entry point / zero-preserving op
        broadcast!(sin, fZ, fX); Z = sparse(fZ)
        broadcast!(sin, Z, X); Z = sparse(fZ) # warmup for @allocated
        @test (@allocated broadcast!(sin, Z, X)) == 0
        @test broadcast!(sin, Z, X) == sparse(broadcast!(sin, fZ, fX))
        # --> test broadcast! entry point / not-zero-preserving op
        broadcast!(cos, fZ, fX); Z = sparse(fZ)
        broadcast!(cos, Z, X); Z = sparse(fZ) # warmup for @allocated
        @test (@allocated broadcast!(cos, Z, X)) == 0
        @test broadcast!(cos, Z, X) == sparse(broadcast!(cos, fZ, fX))
        # --> test shape checks for broadcast! entry point
        # TODO strengthen this test, avoiding dependence on checking whether
        # broadcast_indices throws to determine whether sparse broadcast should throw
        try
            Base.Broadcast.check_broadcast_indices(indices(Z), spzeros((shapeX .- 1)...))
        catch
            @test_throws DimensionMismatch broadcast!(sin, Z, spzeros((shapeX .- 1)...))
        end
    end
    # --> test with vector destination (V/fV)
    fV = Array(first(vecs))
    for Xo in vecs # vector target
        X = SparseVector{Float32,Int32}(Xo)
        shapeX, fX = size(X), Array(X)
        # --> test broadcast! entry point / zero-preserving op
        broadcast!(sin, fV, fX); V = sparse(fV)
        broadcast!(sin, V, X); V = sparse(fV) # warmup for @allocated
        @test (@allocated broadcast!(sin, V, X)) == 0
        @test broadcast!(sin, V, X) == sparse(broadcast!(sin, fV, fX))
        # --> test broadcast! entry point / not-zero-preserving
        broadcast!(cos, fV, fX); V = sparse(fV)
        broadcast!(cos, V, X); V = sparse(fV) # warmup for @allocated
        @test (@allocated broadcast!(cos, V, X)) == 0
        @test broadcast!(cos, V, X) == sparse(broadcast!(cos, fV, fX))
        # --> test shape checks for broadcast! entry point
        # TODO strengthen this test, avoiding dependence on checking whether
        # broadcast_indices throws to determine whether sparse broadcast should throw
        try
            Base.Broadcast.check_broadcast_indices(indices(V), spzeros((shapeX .- 1)...))
        catch
            @test_throws DimensionMismatch broadcast!(sin, V, spzeros((shapeX .- 1)...))
        end
    end
    # Tests specific to #19895, i.e. for broadcast!(identity, C, A) specializations
    Z = copy(first(mats)); fZ = Array(Z)
    V = copy(first(vecs)); fV = Array(V)
    for X in (mats..., vecs...)
        @test broadcast!(identity, Z, X) == sparse(broadcast!(identity, fZ, Array(X)))
        X isa SparseVector && @test broadcast!(identity, V, X) == sparse(broadcast!(identity, fV, Array(X)))
    end
end

@testset "broadcast[!] implementation specialized for pairs of (input) sparse vectors/matrices" begin
    N, M, p = 10, 12, 0.3
    f(x, y) = x + y + 1
    mats = (sprand(N, M, p), sprand(N, 1, p), sprand(1, M, p), sprand(1, 1, 1.0), spzeros(1, 1))
    vecs = (sprand(N, p), sprand(1, 1.0), spzeros(1))
    tens = (mats..., vecs...)
    fZ = Array(first(mats))
    for Xo in tens
        X = ndims(Xo) == 1 ? SparseVector{Float32,Int32}(Xo) : SparseMatrixCSC{Float32,Int32}(Xo)
        # use different types to check internal type stability via allocation tests below
        shapeX, fX = size(X), Array(X)
        for Y in tens
            fY = Array(Y)
            # --> test broadcast entry point
            @test broadcast(+, X, Y) == sparse(broadcast(+, fX, fY))
            @test broadcast(*, X, Y) == sparse(broadcast(*, fX, fY))
            @test broadcast(f, X, Y) == sparse(broadcast(f, fX, fY))
            # TODO strengthen this test, avoiding dependence on checking whether
            # broadcast_indices throws to determine whether sparse broadcast should throw
            try
                Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y)
            catch
                @test_throws DimensionMismatch broadcast(+, spzeros((shapeX .- 1)...), Y)
            end
            # --> test broadcast! entry point / +-like zero-preserving op
            broadcast!(+, fZ, fX, fY); Z = sparse(fZ)
            broadcast!(+, Z, X, Y); Z = sparse(fZ) # warmup for @allocated
            @test (@allocated broadcast!(+, Z, X, Y)) == 0
            @test broadcast!(+, Z, X, Y) == sparse(broadcast!(+, fZ, fX, fY))
            # --> test broadcast! entry point / *-like zero-preserving op
            broadcast!(*, fZ, fX, fY); Z = sparse(fZ)
            broadcast!(*, Z, X, Y); Z = sparse(fZ) # warmup for @allocated
            @test (@allocated broadcast!(*, Z, X, Y)) == 0
            @test broadcast!(*, Z, X, Y) == sparse(broadcast!(*, fZ, fX, fY))
            # --> test broadcast! entry point / not zero-preserving op
            broadcast!(f, fZ, fX, fY); Z = sparse(fZ)
            broadcast!(f, Z, X, Y); Z = sparse(fZ) # warmup for @allocated
            @test (@allocated broadcast!(f, Z, X, Y)) == 0
            @test broadcast!(f, Z, X, Y) == sparse(broadcast!(f, fZ, fX, fY))
            # --> test shape checks for both broadcast and broadcast! entry points
            # TODO strengthen this test, avoiding dependence on checking whether
            # broadcast_indices throws to determine whether sparse broadcast should throw
            try
                Base.Broadcast.check_broadcast_indices(indices(Z), spzeros((shapeX .- 1)...), Y)
            catch
                @test_throws DimensionMismatch broadcast!(f, Z, spzeros((shapeX .- 1)...), Y)
            end
        end
    end
end

@testset "broadcast[!] implementation capable of handling >2 (input) sparse vectors/matrices" begin
    N, M, p = 10, 12, 0.3
    f(x, y, z) = x + y + z + 1
    mats = (sprand(N, M, p), sprand(N, 1, p), sprand(1, M, p), sprand(1, 1, 1.0), spzeros(1, 1))
    vecs = (sprand(N, p), sprand(1, 1.0), spzeros(1))
    tens = (mats..., vecs...)
    for Xo in tens
        X = ndims(Xo) == 1 ? SparseVector{Float32,Int32}(Xo) : SparseMatrixCSC{Float32,Int32}(Xo)
        # use different types to check internal type stability via allocation tests below
        shapeX, fX = size(X), Array(X)
        for Y in tens, Z in tens
            fY, fZ = Array(Y), Array(Z)
            # --> test broadcast entry point
            @test broadcast(+, X, Y, Z) == sparse(broadcast(+, fX, fY, fZ))
            @test broadcast(*, X, Y, Z) == sparse(broadcast(*, fX, fY, fZ))
            @test broadcast(f, X, Y, Z) == sparse(broadcast(f, fX, fY, fZ))
            # TODO strengthen this test, avoiding dependence on checking whether
            # broadcast_indices throws to determine whether sparse broadcast should throw
            try
                Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y, Z)
            catch
                @test_throws DimensionMismatch broadcast(+, spzeros((shapeX .- 1)...), Y, Z)
            end
            # --> test broadcast! entry point / +-like zero-preserving op
            fQ = broadcast(+, fX, fY, fZ); Q = sparse(fQ)
            broadcast!(+, Q, X, Y, Z); Q = sparse(fQ) # warmup for @allocated
            @test (@allocated broadcast!(+, Q, X, Y, Z)) == 0
            @test broadcast!(+, Q, X, Y, Z) == sparse(broadcast!(+, fQ, fX, fY, fZ))
            # --> test broadcast! entry point / *-like zero-preserving op
            fQ = broadcast(*, fX, fY, fZ); Q = sparse(fQ)
            broadcast!(*, Q, X, Y, Z); Q = sparse(fQ) # warmup for @allocated
            @test (@allocated broadcast!(*, Q, X, Y, Z)) == 0
            @test broadcast!(*, Q, X, Y, Z) == sparse(broadcast!(*, fQ, fX, fY, fZ))
            # --> test broadcast! entry point / not zero-preserving op
            fQ = broadcast(f, fX, fY, fZ); Q = sparse(fQ)
            broadcast!(f, Q, X, Y, Z); Q = sparse(fQ) # warmup for @allocated
            @test_broken (@allocated broadcast!(f, Q, X, Y, Z)) == 0
            # the preceding test allocates 16 bytes in the entry point for broadcast!, but
            # none of the earlier tests of the same code path allocate. no allocation shows
            # up with --track-allocation=user. allocation shows up on the first line of the
            # entry point for broadcast! with --track-allocation=all, but that first line
            # almost certainly should not allocate. so not certain what's going on.
            # additional info: occurs for broadcast!(f, Z, X) for Z and X of different
            # shape, but not for Z and X of the same shape.
            @test broadcast!(f, Q, X, Y, Z) == sparse(broadcast!(f, fQ, fX, fY, fZ))
            # --> test shape checks for both broadcast and broadcast! entry points
            # TODO strengthen this test, avoiding dependence on checking whether
            # broadcast_indices throws to determine whether sparse broadcast should throw
            try
                Base.Broadcast.check_broadcast_indices(indices(Q), spzeros((shapeX .- 1)...), Y, Z)
            catch
                @test_throws DimensionMismatch broadcast!(f, Q, spzeros((shapeX .- 1)...), Y, Z)
            end
        end
    end
end

@testset "sparse map/broadcast with result eltype not a concrete subtype of Number (#19561/#19589)" begin
    N = 4
    A, fA = speye(N), eye(N)
    B, fB = spzeros(1, N), zeros(1, N)
    intorfloat_zeropres(xs...) = all(iszero, xs) ? zero(Float64) : Int(1)
    stringorfloat_zeropres(xs...) = all(iszero, xs) ? zero(Float64) : "hello"
    intorfloat_notzeropres(xs...) = all(iszero, xs) ? Int(1) : zero(Float64)
    stringorfloat_notzeropres(xs...) = all(iszero, xs) ? "hello" : zero(Float64)
    for fn in (intorfloat_zeropres, intorfloat_notzeropres,
                stringorfloat_zeropres, stringorfloat_notzeropres)
        @test map(fn, A) == sparse(map(fn, fA))
        @test broadcast(fn, A) == sparse(broadcast(fn, fA))
        @test broadcast(fn, A, B) == sparse(broadcast(fn, fA, fB))
        @test broadcast(fn, B, A) == sparse(broadcast(fn, fB, fA))
    end
    for fn in (intorfloat_zeropres, stringorfloat_zeropres)
        @test broadcast(fn, A, B, A) == sparse(broadcast(fn, fA, fB, fA))
    end
end

@testset "broadcast[!] over combinations of scalars and sparse vectors/matrices" begin
    N, M, p = 10, 12, 0.5
    elT = Float64
    s = Float32(2.0)
    V = sprand(elT, N, p)
    A = sprand(elT, N, M, p)
    fV, fA = Array(V), Array(A)
    # test combinations involving one to three scalars and one to five sparse vectors/matrices
    spargseq, dargseq = Iterators.cycle((A, V)), Iterators.cycle((fA, fV))
    for nargs in 1:5 # number of tensor arguments
        nargsl = cld(nargs, 2) # number in "left half" of tensor arguments
        nargsr = fld(nargs, 2) # number in "right half" of tensor arguments
        spargsl = tuple(Iterators.take(spargseq, nargsl)...) # "left half" of tensor args
        spargsr = tuple(Iterators.take(spargseq, nargsr)...) # "right half" of tensor args
        dargsl = tuple(Iterators.take(dargseq, nargsl)...) # "left half" of tensor args, densified
        dargsr = tuple(Iterators.take(dargseq, nargsr)...) # "right half" of tensor args, densified
        for (sparseargs, denseargs) in ( # argument combinations including scalars
                # a few combinations involving one scalar
                ((s, spargsl..., spargsr...), (s, dargsl..., dargsr...)),
                ((spargsl..., s, spargsr...), (dargsl..., s, dargsr...)),
                ((spargsl..., spargsr..., s), (dargsl..., dargsr..., s)),
                # a few combinations involving two scalars
                ((s, spargsl..., s, spargsr...), (s, dargsl..., s, dargsr...)),
                ((s, spargsl..., spargsr..., s), (s, dargsl..., dargsr..., s)),
                ((spargsl..., s, spargsr..., s), (dargsl..., s, dargsr..., s)),
                ((s, s, spargsl..., spargsr...), (s, s, dargsl..., dargsr...)),
                ((spargsl..., s, s, spargsr...), (dargsl..., s, s, dargsr...)),
                ((spargsl..., spargsr..., s, s), (dargsl..., dargsr..., s, s)),
                # a few combinations involving three scalars
                ((s, spargsl..., s, spargsr..., s), (s, dargsl..., s, dargsr..., s)),
                ((s, spargsl..., s, s, spargsr...), (s, dargsl..., s, s, dargsr...)),
                ((spargsl..., s, s, spargsr..., s), (dargsl..., s, s, dargsr..., s)),
                ((spargsl..., s, s, s, spargsr...), (dargsl..., s, s, s, dargsr...)), )
            # test broadcast entry point
            @test broadcast(*, sparseargs...) == sparse(broadcast(*, denseargs...))
            @test isa(@inferred(broadcast(*, sparseargs...)), SparseMatrixCSC{elT})
            # test broadcast! entry point
            fX = broadcast(*, sparseargs...); X = sparse(fX)
            @test broadcast!(*, X, sparseargs...) == sparse(broadcast!(*, fX, denseargs...))
            @test isa(@inferred(broadcast!(*, X, sparseargs...)), SparseMatrixCSC{elT})
            X = sparse(fX) # reset / warmup for @allocated test
            @test_broken (@allocated broadcast!(*, X, sparseargs...)) == 0
            # This test (and the analog below) fails for three reasons:
            # (1) In all cases, generating the closures that capture the scalar arguments
            #   results in allocation, not sure why.
            # (2) In some cases, though _broadcast_eltype (which wraps _return_type)
            #   consistently provides the correct result eltype when passed the closure
            #   that incorporates the scalar arguments to broadcast (and, with #19667,
            #   is inferable, so the overall return type from broadcast is inferred),
            #   in some cases inference seems unable to determine the return type of
            #   direct calls to that closure. This issue causes variables in both the
            #   broadcast[!] entry points (fofzeros = f(_zeros_eltypes(args...)...)) and
            #   the driver routines (Cx in _map_zeropres! and _broadcast_zeropres!) to have
            #   inferred type Any, resulting in allocation and lackluster performance.
            # (3) The sparseargs... splat in the call above allocates a bit, but of course
            #   that issue is negligible and perhaps could be accounted for in the test.
        end
    end
    # test combinations at the limit of inference (eight arguments net)
    for (sparseargs, denseargs) in (
            ((s, s, s, A, s, s, s, s), (s, s, s, fA, s, s, s, s)), # seven scalars, one sparse matrix
            ((s, s, V, s, s, A, s, s), (s, s, fV, s, s, fA, s, s)), # six scalars, two sparse vectors/matrices
            ((s, s, V, s, A, s, V, s), (s, s, fV, s, fA, s, fV, s)), # five scalars, three sparse vectors/matrices
            ((s, V, s, A, s, V, s, A), (s, fV, s, fA, s, fV, s, fA)), # four scalars, four sparse vectors/matrices
            ((s, V, A, s, V, A, s, A), (s, fV, fA, s, fV, fA, s, fA)), # three scalars, five sparse vectors/matrices
            ((V, A, V, s, A, V, A, s), (fV, fA, fV, s, fA, fV, fA, s)), # two scalars, six sparse vectors/matrices
            ((V, A, V, A, s, V, A, V), (fV, fA, fV, fA, s, fV, fA, fV)) ) # one scalar, seven sparse vectors/matrices
        # test broadcast entry point
        @test broadcast(*, sparseargs...) == sparse(broadcast(*, denseargs...))
        @test isa(@inferred(broadcast(*, sparseargs...)), SparseMatrixCSC{elT})
        # test broadcast! entry point
        fX = broadcast(*, sparseargs...); X = sparse(fX)
        @test broadcast!(*, X, sparseargs...) == sparse(broadcast!(*, fX, denseargs...))
        @test isa(@inferred(broadcast!(*, X, sparseargs...)), SparseMatrixCSC{elT})
        X = sparse(fX) # reset / warmup for @allocated test
        @test_broken (@allocated broadcast!(*, X, sparseargs...)) == 0
        # please see the note a few lines above re. this @test_broken
    end
end

@testset "broadcast[!] over combinations of scalars, sparse arrays, structured matrices, and dense vectors/matrices" begin
    N, p = 10, 0.4
    s = rand()
    V = sprand(N, p)
    A = sprand(N, N, p)
    Z = copy(A)
    sparsearrays = (V, A)
    fV, fA = map(Array, sparsearrays)
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), true)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    S = SymTridiagonal(rand(N), rand(N - 1))
    structuredarrays = (D, B, T, S)
    fstructuredarrays = map(Array, structuredarrays)
    for (X, fX) in zip(structuredarrays, fstructuredarrays)
        @test (Q = broadcast(sin, X); Q isa SparseMatrixCSC && Q == sparse(broadcast(sin, fX)))
        @test broadcast!(sin, Z, X) == sparse(broadcast(sin, fX))
        @test (Q = broadcast(cos, X); Q isa SparseMatrixCSC && Q == sparse(broadcast(cos, fX)))
        @test broadcast!(cos, Z, X) == sparse(broadcast(cos, fX))
        @test (Q = broadcast(*, s, X); Q isa SparseMatrixCSC && Q == sparse(broadcast(*, s, fX)))
        @test broadcast!(*, Z, s, X) == sparse(broadcast(*, s, fX))
        @test (Q = broadcast(+, V, A, X); Q isa SparseMatrixCSC && Q == sparse(broadcast(+, fV, fA, fX)))
        @test broadcast!(+, Z, V, A, X) == sparse(broadcast(+, fV, fA, fX))
        @test (Q = broadcast(*, s, V, A, X); Q isa SparseMatrixCSC && Q == sparse(broadcast(*, s, fV, fA, fX)))
        @test broadcast!(*, Z, s, V, A, X) == sparse(broadcast(*, s, fV, fA, fX))
        for (Y, fY) in zip(structuredarrays, fstructuredarrays)
            @test (Q = broadcast(+, X, Y); Q isa SparseMatrixCSC && Q == sparse(broadcast(+, fX, fY)))
            @test broadcast!(+, Z, X, Y) == sparse(broadcast(+, fX, fY))
            @test (Q = broadcast(*, X, Y); Q isa SparseMatrixCSC && Q == sparse(broadcast(*, fX, fY)))
            @test broadcast!(*, Z, X, Y) == sparse(broadcast(*, fX, fY))
        end
    end
    C = Array(sprand(N, 0.4))
    M = Array(sprand(N, N, 0.4))
    densearrays = (C, M)
    fD, fB = Array(D), Array(B)
    for X in densearrays
        @test broadcast(+, D, X)::SparseMatrixCSC == sparse(broadcast(+, fD, X))
        @test broadcast!(+, Z, D, X) == sparse(broadcast(+, fD, X))
        @test broadcast(*, s, B, X)::SparseMatrixCSC == sparse(broadcast(*, s, fB, X))
        @test broadcast!(*, Z, s, B, X) == sparse(broadcast(*, s, fB, X))
        @test broadcast(+, V, B, X)::SparseMatrixCSC == sparse(broadcast(+, fV, fB, X))
        @test broadcast!(+, Z, V, B, X) == sparse(broadcast(+, fV, fB, X))
        @test broadcast(+, V, A, X)::SparseMatrixCSC == sparse(broadcast(+, fV, fA, X))
        @test broadcast!(+, Z, V, A, X) == sparse(broadcast(+, fV, fA, X))
        @test broadcast(*, s, V, A, X)::SparseMatrixCSC == sparse(broadcast(*, s, fV, fA, X))
        @test broadcast!(*, Z, s, V, A, X) == sparse(broadcast(*, s, fV, fA, X))
        # Issue #20954 combinations of sparse arrays and RowVectors
        @test broadcast(+, A, X')::SparseMatrixCSC == sparse(broadcast(+, fA, X'))
        @test broadcast(*, V, X')::SparseMatrixCSC == sparse(broadcast(*, fV, X'))
    end
end

@testset "broadcast! where the destination is a structured matrix" begin
    # Where broadcast!'s destination is a structured matrix, broadcast! should fall back
    # to the generic AbstractArray broadcast! code (at least for now).
    N, p = 5, 0.4
    A = sprand(N, N, p)
    sA = A + transpose(A)
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), true)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D))
    @test broadcast!(sin, copy(B), B) == Bidiagonal(sin.(B), true)
    @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T))
    @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A))
    @test broadcast!(*, copy(B), B, A) == Bidiagonal(broadcast(*, B, A), true)
    @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A))
    # SymTridiagonal (and similar symmetric matrix types) do not support setindex!
    # off the diagonal, and so cannot serve as a destination for broadcast!
end

@testset "map[!] over combinations of sparse and structured matrices" begin
    N, p = 10, 0.4
    A = sprand(N, N, p)
    Z, fA = copy(A), Array(A)
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), true)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    S = SymTridiagonal(rand(N), rand(N - 1))
    structuredarrays = (D, B, T, S)
    fstructuredarrays = map(Array, structuredarrays)
    for (X, fX) in zip(structuredarrays, fstructuredarrays)
        @test (Q = map(sin, X); Q isa SparseMatrixCSC && Q == sparse(map(sin, fX)))
        @test map!(sin, Z, X) == sparse(map(sin, fX))
        @test (Q = map(cos, X); Q isa SparseMatrixCSC && Q == sparse(map(cos, fX)))
        @test map!(cos, Z, X) == sparse(map(cos, fX))
        @test (Q = map(+, A, X); Q isa SparseMatrixCSC && Q == sparse(map(+, fA, fX)))
        @test map!(+, Z, A, X) == sparse(map(+, fA, fX))
        for (Y, fY) in zip(structuredarrays, fstructuredarrays)
            @test (Q = map(+, X, Y); Q isa SparseMatrixCSC && Q == sparse(map(+, fX, fY)))
            @test map!(+, Z, X, Y) == sparse(map(+, fX, fY))
            @test (Q = map(*, X, Y); Q isa SparseMatrixCSC && Q == sparse(map(*, fX, fY)))
            @test map!(*, Z, X, Y) == sparse(map(*, fX, fY))
            @test (Q = map(+, X, A, Y); Q isa SparseMatrixCSC && Q == sparse(map(+, fX, fA, fY)))
            @test map!(+, Z, X, A, Y) == sparse(map(+, fX, fA, fY))
        end
    end
end

# Older tests of sparse broadcast, now largely covered by the tests above
@testset "assorted tests of sparse broadcast over two input arguments" begin
    N, p = 10, 0.3
    A, B, CF = sprand(N, N, p), sprand(N, N, p), rand(N, N)
    AF, BF, C = Array(A), Array(B), sparse(CF)

    @test A .* B == AF .* BF
    @test A[1,:] .* B == AF[1,:] .* BF
    @test A[:,1] .* B == AF[:,1] .* BF
    @test A .* B[1,:] == AF .*  BF[1,:]
    @test A .* B[:,1] == AF .*  BF[:,1]

    @test A .* B == AF .* BF
    @test A[1,:] .* BF == AF[1,:] .* BF
    @test A[:,1] .* BF == AF[:,1] .* BF
    @test A .* BF[1,:] == AF .*  BF[1,:]
    @test A .* BF[:,1] == AF .*  BF[:,1]

    @test A .* B == AF .* BF
    @test AF[1,:] .* B == AF[1,:] .* BF
    @test AF[:,1] .* B == AF[:,1] .* BF
    @test AF .* B[1,:] == AF .*  BF[1,:]
    @test AF .* B[:,1] == AF .*  BF[:,1]

    @test A .* B == AF .* BF
    @test A[1,:] .* B == AF[1,:] .* BF
    @test A[:,1] .* B == AF[:,1] .* BF
    @test A .* B[1,:] == AF .*  BF[1,:]
    @test A .* B[:,1] == AF .*  BF[:,1]

    @test A .* 3 == AF .* 3
    @test 3 .* A == 3 .* AF
    @test A[1,:] .* 3 == AF[1,:] .* 3
    @test A[:,1] .* 3 == AF[:,1] .* 3

    @test A .- 3 == AF .- 3
    @test 3 .- A == 3 .- AF
    @test A .- B == AF .- BF
    @test A - AF == zeros(AF)
    @test AF - A == zeros(AF)
    @test A[1,:] .- B == AF[1,:] .- BF
    @test A[:,1] .- B == AF[:,1] .- BF
    @test A .- B[1,:] == AF .-  BF[1,:]
    @test A .- B[:,1] == AF .-  BF[:,1]

    @test A .+ 3 == AF .+ 3
    @test 3 .+ A == 3 .+ AF
    @test A .+ B == AF .+ BF
    @test A + AF == AF + A
    @test (A .< B) == (AF .< BF)
    @test (A .!= B) == (AF .!= BF)

    @test A ./ 3 == AF ./ 3
    @test A .\ 3 == AF .\ 3
    @test 3 ./ A == 3 ./ AF
    @test 3 .\ A == 3 .\ AF
    @test A .\ C == AF .\ CF
    @test A ./ C == AF ./ CF
    @test A ./ CF[:,1] == AF ./ CF[:,1]
    @test A .\ CF[:,1] == AF .\ CF[:,1]
    @test BF ./ C == BF ./ CF
    @test BF .\ C == BF .\ CF

    @test A .^ 3 == AF .^ 3
    @test 3 .^ A == 3 .^ AF
    @test A .^ BF[:,1] == AF .^ BF[:,1]
    @test BF[:,1] .^ A == BF[:,1] .^ AF

    @test spzeros(0,0)  + spzeros(0,0) == zeros(0,0)
    @test spzeros(0,0)  * spzeros(0,0) == zeros(0,0)
    @test spzeros(1,0) .+ spzeros(2,1) == zeros(2,0)
    @test spzeros(1,0) .* spzeros(2,1) == zeros(2,0)
    @test spzeros(1,2) .+ spzeros(0,1) == zeros(0,2)
    @test spzeros(1,2) .* spzeros(0,1) == zeros(0,2)
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using  IterativeSolvers, Roots # , KrylovKit
export solve_ω, _solve_Δω², solve_k, solve_n, ng, k_guess, solve_nω, solve_ω², make_εₛ⁻¹, make_MG, make_MD

"""
################################################################################
#																			   #
#	Routines to shield expensive initialization calculations from memory-	   #
#						intensive reverse-mode AD				  			   #
#																			   #
################################################################################
"""

k_guess(ω,ε⁻¹::Array{Float64,5}) = ( kg = Zygote.@ignore ( first(ω) * sqrt(1/minimum([minimum(ε⁻¹[a,a,:,:,:]) for a=1:3])) ); kg  )
k_guess(ω,shapes::Vector{<:Shape}) = ( kg = Zygote.@ignore ( first(ω) * √εₘₐₓ(shapes) ); kg  )
make_MG(Δx,Δy,Δz,Nx,Ny,Nz) = (g = Zygote.@ignore (MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)); g)::MaxwellGrid
make_MD(k,g::MaxwellGrid) = (ds = Zygote.@ignore (MaxwellData(k,g)); ds)::MaxwellData
make_KDTree(shapes::Vector{<:Shape}) = (tree = Zygote.@ignore (KDTree(shapes)); tree)::KDTree

function make_εₛ⁻¹(shapes::Vector{<:Shape},g::MaxwellGrid)::Array{Float64,5}
    tree = make_KDTree(shapes)
    eibuf = Zygote.Buffer(Array{Float64}(undef),3,3,g.Nx,g.Ny,1)
    for i=1:g.Nx,j=1:g.Ny,kk=1:g.Nz
		# eps = εₛ(shapes,Zygote.dropgrad(tree),Zygote.dropgrad(g.x[i]),Zygote.dropgrad(g.y[j]),Zygote.dropgrad(g.δx),Zygote.dropgrad(g.δy))
		eps = εₛ(shapes,tree,g.x[i],g.y[j],g.δx,g.δy)
		epsi = inv(eps) # inv( (eps' + eps) / 2) # Hermitian(inv(eps))  # inv(Hermitian(eps)) #   # inv(eps)
        eibuf[:,:,i,j,kk] = epsi #(epsi' + epsi) / 2
    end
    return real(copy(eibuf))
end

"""
################################################################################
#																			   #
#						solve_ω² methods: (ε⁻¹, k) --> (H, ω²)				   #
#																			   #
################################################################################
"""

function solve_ω²(kz::T,ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8) where T<:Real
	# Δk = k - ds.k
	ds.k = kz
	ds.kpg_mag, ds.mn = calc_kpg(kz,ds.Δx,ds.Δy,ds.Δz,ds.Nx,ds.Ny,ds.Nz)
    # res = IterativeSolvers.lobpcg(M̂(ε⁻¹,ds),false,neigs;P=P̂(ε⁻¹,ds),maxiter,tol)
    res = IterativeSolvers.lobpcg(M̂!(ε⁻¹,ds),false,ds.H⃗;P=P̂!(ε⁻¹,ds),maxiter,tol)
    H =  res.X #[:,eigind]                       # eigenmode wavefn. magnetic fields in transverse pol. basis
    ds.ω² =  real(res.λ[eigind])                     # eigenmode temporal freq.,  neff = kz / ω, kz = k[3]
	ds.H⃗ .= H
	ds.ω = ( ds.ω² > 0. ? sqrt(ds.ω²) : 0. )
    # ds.ω²ₖ = 2 * H_Mₖ_H(Ha,ε⁻¹,kpg_mn,kpg_mag,ds.𝓕,ds.𝓕⁻¹) # = 2ω*ωₖ; ωₖ = ∂ω/∂kz = group velocity = c / ng; c = 1 here
    return (H, ds.ω^2) #(H, ds.ω²) #, ωₖ
end
# @btime: solve_ω²(1.5,$ε⁻¹_mpb;ds=$ds)
# 536.372 ms (17591 allocations: 125.75 MiB)

# function solve_ω²(k::Array{<:Real},ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	outs = [solve_ω²(kk,ε⁻¹,ds;neigs,eigind,maxiter,tol) for kk in k]
#     ( [o[1] for o in outs], [o[2] for o in outs] )
# end

function solve_ω²(kz,ε⁻¹::Array{Float64,5},g::MaxwellGrid;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    solve_ω²(kz,ε⁻¹,make_MD(first(kz),g);neigs,eigind,maxiter,tol)
end
# @btime:
# 498.442 ms (13823 allocations: 100.19 MiB)

function solve_ω²(kz::T,ε⁻¹::Array{T,5},Δx::T,Δy::T,Δz::T;neigs=1,eigind=1,maxiter=3000,tol=1e-8) where T<:Real
    solve_ω²(kz,ε⁻¹,make_MG(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...);neigs,eigind,maxiter,tol)
end

# function solve_k(ω::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},g::MaxwellGrid;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
#     solve_k(ω,shapes,make_MD(kguess,g)::MaxwellData;neigs,eigind,maxiter,tol)
# end

# function solve_k(ω::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	g = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
#     solve_k(ω,shapes,g;kguess,neigs,eigind,maxiter,tol)
# end

function solve_ω²(kz,shapes::Vector{<:GeometryPrimitives.Shape},g::MaxwellGrid;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	ds::MaxwellData = make_MD(kz,g)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	solve_ω²(kz,ε⁻¹,ds;neigs,eigind,maxiter,tol)
end

function solve_ω²(kz,shapes::Vector{<:GeometryPrimitives.Shape},Δx,Δy,Δz,Nx,Ny,Nz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ds::MaxwellData = make_MD(kz,g)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	solve_ω²(kz,ε⁻¹,ds;neigs,eigind,maxiter,tol)
end

# function solve_k(ω::Number,shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
# 	ds::MaxwellData = make_MD(kguess,g)
# 	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
#     # solve_k(ω,ε⁻¹,ds;neigs,eigind,maxiter,tol)
# 	kz = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), ds.k, Roots.Newton())
# 	return ( copy(ds.H⃗), kz )
# end

# function solve_ω²(k::Vector{<:Number},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
# 	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
#     outs = [solve_ω²(kk,ε⁻¹,make_MD(kk,g)::MaxwellData;neigs,eigind,maxiter,tol) for kk in k]
# 	return ( [o[1] for o in outs], [o[2] for o in outs] ) #( copy(ds.H⃗), kz )
# end

"""
################################################################################
#																			   #
#						solve_ω methods: (ε⁻¹, k) --> (H, ω)				   #
#																			   #
################################################################################
"""

function solve_ω(k::T,ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8) where T<:Real
	# Δk = k - ds.k
	ds.k = k
	ds.kpg_mag, ds.mn = calc_kpg(k,ds.Δx,ds.Δy,ds.Δz,ds.Nx,ds.Ny,ds.Nz)
    # res = IterativeSolvers.lobpcg(M̂(ε⁻¹,ds),false,neigs;P=P̂(ε⁻¹,ds),maxiter,tol)
    res = IterativeSolvers.lobpcg(M̂!(ε⁻¹,ds),false,ds.H⃗;P=P̂!(ε⁻¹,ds),maxiter,tol)
    H =  res.X #[:,eigind]                       # eigenmode wavefn. magnetic fields in transverse pol. basis
    ω =  √(real(res.λ[eigind]))                     # eigenmode temporal freq.,  neff = kz / ω, kz = k[3]
	ds.H⃗ .= H
    ds.ω² = ω^2; ds.ω = ω;
    # ds.ω²ₖ = 2 * H_Mₖ_H(Ha,ε⁻¹,kpg_mn,kpg_mag,ds.𝓕,ds.𝓕⁻¹) # = 2ω*ωₖ; ωₖ = ∂ω/∂kz = group velocity = c / ng; c = 1 here
    return H, ω #, ωₖ
end
# @btime: solve_ω(1.5,$ε⁻¹_mpb;ds=$ds)
# 536.372 ms (17591 allocations: 125.75 MiB)

function solve_ω(k::Array{<:Real},ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	outs = [solve_ω(kk,ε⁻¹,ds;neigs,eigind,maxiter,tol) for kk in k]
    ( [o[1] for o in outs], [o[2] for o in outs] )
end

function solve_ω(k,ε⁻¹::Array{Float64,5},g::MaxwellGrid;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    solve_ω(k,ε⁻¹,MaxwellData(first(k),g);neigs,eigind,maxiter,tol)
end
# @btime:
# 498.442 ms (13823 allocations: 100.19 MiB)

function solve_ω(k,ε⁻¹::AbstractArray,Δx,Δy,Δz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    solve_ω(k,ε⁻¹,MaxwellGrid(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...);neigs,eigind,maxiter,tol)
end

# function solve_ω(k,shapes::Vector{<:Shape},Δx,Δy,Δz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
#     solve_ω(k,ε⁻¹,MaxwellGrid(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...);neigs,eigind,maxiter,tol)
# end
#
# function solve_ω(k::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
#     solve_k(ω,make_εₛ⁻¹(shapes,ds.grid)::Array{Float64,5},ds;neigs,eigind,maxiter,tol)
# end

"""
################################################################################
#																			   #
#						solve_k methods: (ε⁻¹, ω) --> (H, k)				   #
#																			   #
################################################################################
"""

"""
modified solve_ω version for Newton solver, which wants (x -> f(x), f(x)/f'(x)) as input to solve f(x) = 0
"""
function _solve_Δω²(k,ωₜ,ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    ds.k = k
	ds.kpg_mag, ds.mn = calc_kpg(k,ds.Δx,ds.Δy,ds.Δz,ds.Nx,ds.Ny,ds.Nz)
    res = IterativeSolvers.lobpcg(M̂!(ε⁻¹,ds),false,ds.H⃗;P=P̂!(ε⁻¹,ds),maxiter,tol)
    ds.H⃗ .=  res.X #[:,eigind]                      # eigenmode wavefn. magnetic fields in transverse pol. basis
    ds.ω² =  (real(res.λ[eigind]))                # eigenmode temporal freq.,  neff = kz / ωₖ, kz = k[3]
    Δω² = ds.ω² - ωₜ^2
    # ω²ₖ =   2 * real( ( H[:,eigind]' * M̂ₖ(ε⁻¹,ds) * H[:,eigind] )[1])  # = 2ω*ωₖ; ωₖ = ∂ω/∂kz = group velocity = c / ng; c = 1 here
	# Ha = reshape(H,(2,size(ε⁻¹)[end-2:end]...))
	ds.ω²ₖ = 2 * H_Mₖ_H(ds.H⃗,ε⁻¹,ds.kpg_mag,ds.mn) #,ds.𝓕,ds.𝓕⁻¹) # = 2ω*ωₖ; ωₖ = ∂ω/∂kz = group velocity = c / ng; c = 1 here
    return Δω² , Δω² / ds.ω²ₖ
end

# function solve_k(ω,ε⁻¹;Δx=6.0,Δy=4.0,Δz=1.0,k_guess=ω*sqrt(1/minimum([minimum(ε⁻¹[a,a,:,:,:]) for a=1:3])),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
function solve_k(ω::Number,ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    kz = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), ds.k, Roots.Newton())
    return ( copy(ds.H⃗), kz ) # maybe copy(ds.H⃗) instead?
end

function solve_k(ω::Vector{<:Number},ε⁻¹::Array{Float64,5},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    outs = [solve_k(om,ε⁻¹,ds;neigs,eigind,maxiter,tol) for om in ω]
    ( [o[1] for o in outs], [o[2] for o in outs] )
end

function solve_k(ω::Number,ε⁻¹::Array{Float64,5},g::MaxwellGrid;kguess=k_guess(ω,ε⁻¹),neigs=1,eigind=1,maxiter=1000,tol=1e-6)
	solve_k(ω,ε⁻¹,make_MD(kguess,g);neigs,eigind,maxiter,tol)
end

function solve_k(ω::Number,ε⁻¹::Array{Float64,5},Δx,Δy,Δz;kguess=k_guess(ω,ε⁻¹),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g = make_MG(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...) #MaxwellGrid(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...)
	ds = make_MD(kguess,g) 			# MaxwellData(kguess,g)
    solve_k(ω,ε⁻¹,ds;neigs,eigind,maxiter,tol)
end

function solve_k(ω::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    solve_k(ω,make_εₛ⁻¹(shapes,ds.grid)::Array{Float64,5},ds;neigs,eigind,maxiter,tol)
end

function solve_k(ω::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},g::MaxwellGrid;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    solve_k(ω,shapes,make_MD(kguess,g)::MaxwellData;neigs,eigind,maxiter,tol)
end

# function solve_k(ω::Union{Number,Vector{<:Number}},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	g = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
#     solve_k(ω,shapes,g;kguess,neigs,eigind,maxiter,tol)
# end

function solve_k(ω::Number,shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ds::MaxwellData = make_MD(kguess,g)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
    # solve_k(ω,ε⁻¹,ds;neigs,eigind,maxiter,tol)
	kz = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), ds.k, Roots.Newton())
	return ( copy(ds.H⃗), kz )
end

function solve_k(ω::Vector{<:Number},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ds::MaxwellData = make_MD(kguess,g)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
    outs = [solve_k(om,ε⁻¹,ds;neigs,eigind,maxiter,tol) for om in ω]
	return ( [o[1] for o in outs], [o[2] for o in outs] ) #( copy(ds.H⃗), kz )
end

"""
################################################################################
#																			   #
#						solve_n methods: (ε⁻¹, ω) --> (n, ng)				   #
#																			   #
################################################################################
"""

function solve_n(ω::Number,ε⁻¹::AbstractArray,ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	k = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), ds.k, Roots.Newton())
	( k / ω , 2ω / ds.ω²ₖ ) # = ( n , ng )
end

function solve_n(ω::Array{<:Real,1},ε⁻¹::AbstractArray,ds::MaxwellData;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
    outs = [solve_n(om,ε⁻¹,ds;neigs,eigind,maxiter,tol) for om in ω]
    ( [o[1] for o in outs], [o[2] for o in outs] )
end

function solve_n(ω,ε⁻¹::Array{Float64,5},g::MaxwellGrid;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	k_guess = first(ω) * sqrt(1/minimum([minimum(ε⁻¹[a,a,:,:,:]) for a=1:3]))
	solve_n(ω,ε⁻¹,MaxwellData(k_guess,g);neigs,eigind,maxiter,tol)
end
# @btime:
# 498.442 ms (13823 allocations: 100.19 MiB)

# function solve_n(ω,ε⁻¹::AbstractArray,Δx,Δy,Δz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
#     solve_n(ω,ε⁻¹,MaxwellGrid(Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...);neigs,eigind,maxiter,tol)
# end

function solve_n(ω::Array{<:Real},ε⁻¹::Array{<:Real,5},Δx::T,Δy::T,Δz::T;eigind=1,maxiter=3000,tol=1e-8) where T<:Real
	H,k = solve_k(ω, ε⁻¹,Δx,Δy,Δz;eigind,maxiter,tol)
	( k ./ ω, [ ω[i] / H_Mₖ_H(H[i],ε⁻¹,calc_kpg(k[i],Δx,Δy,Δz,size(ε⁻¹)[end-2:end]...)...) for i=1:length(ω) ] ) # = (n, ng)
end

function solve_n(ω,ε⁻¹,Δx,Δy,Δz;eigind=1,maxiter=3000,tol=1e-8)
	Nx,Ny,Nz = size(ε⁻¹)[end-2:end]
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)
	H,kz = solve_k(ω, ε⁻¹,Δx,Δy,Δz)
	mag, mn = calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	ng = ω / H_Mₖ_H(H,ε⁻¹,mag,mn)
	( kz/ω, ng )
end

# function solve_n(ω::Array{<:Real},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	H,k = solve_k(ω,shapes,Δx,Δy,Δz,Nx,Ny,Nz;kguess,neigs,eigind,maxiter,tol)
# 	( k ./ ω, [ ω[i] / H_Mₖ_H(H[i],ε⁻¹,calc_kpg(k[i],Δx,Δy,Δz,Nx,Ny,Nz)...) for i=1:length(ω) ] ) # = (n, ng)
# end

# function solve_n(ω::Number,shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
# 	H,k = solve_k(ω,shapes,Δx,Δy,Δz,Nx,Ny,Nz;kguess,neigs,eigind,maxiter,tol)
# 	ng = ω / H_Mₖ_H(H,ε⁻¹,calc_kpg(k,Δx,Δy,Δz,Nx,Ny,Nz)...)
# 	( k/ω, ng )
# end

function solve_n(ω::Number,shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	H,kz = solve_k(ω,ε⁻¹,Δx,Δy,Δz;kguess,neigs,eigind,maxiter,tol)
	kpg_mag,kpg_mn = calc_kpg(kz,Zygote.dropgrad(Δx),Zygote.dropgrad(Δy),Zygote.dropgrad(Δz),Zygote.dropgrad(Nx),Zygote.dropgrad(Ny),Zygote.dropgrad(Nz))
	ng = ω / H_Mₖ_H(H,ε⁻¹,kpg_mag,kpg_mn)
	( kz/ω, ng )
end

function solve_n(ω::Array{<:Real},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;kguess=k_guess(ω,shapes),neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	H,k = solve_k(ω,shapes,Δx,Δy,Δz,Nx,Ny,Nz;kguess,neigs,eigind,maxiter,tol)
	( k ./ ω, [ ω[i] / H_Mₖ_H(H[i],ε⁻¹,calc_kpg(k[i],Δx,Δy,Δz,Nx,Ny,Nz)...) for i=1:length(ω) ] ) # = (n, ng)
end

"""
################################################################################
#																			   #
#					solve_nω methods: (ε⁻¹, k) --> (n, ng)					   #
#						(mostly for debugging gradients)					   #
#																			   #
################################################################################
"""

function solve_nω(kz::T,shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;neigs=1,eigind=1,maxiter=3000,tol=1e-8) where T<:Real
	# g::MaxwellGrid = make_MG(Zygote.dropgrad(Δx),Zygote.dropgrad(Δy),Zygote.dropgrad(Δz),Zygote.dropgrad(Nx),Zygote.dropgrad(Ny),Zygote.dropgrad(Nz)) #Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	# ds::MaxwellData = make_MD(kz,g) # MaxwellData(kz,g)
	# kpg_mag,kpg_mn = calc_kpg(kz,Zygote.dropgrad(Δx),Zygote.dropgrad(Δy),Zygote.dropgrad(Δz),Zygote.dropgrad(Nx),Zygote.dropgrad(Ny),Zygote.dropgrad(Nz))
	# mag,mn = calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	mag,mn = calc_kpg(kz,g.g⃗)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	H,ω² = solve_ω²(kz,ε⁻¹,Δx,Δy,Δz;neigs,eigind,maxiter,tol)
	# println("ω² = $ω²")
	ω = sqrt(ω²)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	# ng = -ω / real( dot(Ha, kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,ds.mn), (2:4) ), ε⁻¹), (2:4)),ds.mn,ds.mag) ) )
	ng = ω / H_Mₖ_H(Ha,ε⁻¹,mag,mn)
	# ng = ω / real( dot(H, -vec( kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,mn), (2:4) ), ε⁻¹), (2:4)),mn,mag) ) ) )
	# ng = -ω / real( dot(Ha, kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,Zygote.@showgrad(mn)), (2:4) ), ε⁻¹), (2:4)), Zygote.@showgrad(mn),Zygote.@showgrad(mag)) ) )
	( kz/ω, ng )
end

function solve_nω(kz::T,ε⁻¹::Array{T,5},Δx,Δy,Δz,Nx,Ny,Nz;neigs=1,eigind=1,maxiter=3000,tol=1e-8) where T<:Real
	# g::MaxwellGrid = make_MG(Zygote.dropgrad(Δx),Zygote.dropgrad(Δy),Zygote.dropgrad(Δz),Zygote.dropgrad(Nx),Zygote.dropgrad(Ny),Zygote.dropgrad(Nz)) #Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	# ds::MaxwellData = make_MD(kz,g) # MaxwellData(kz,g)
	# kpg_mag,kpg_mn = calc_kpg(kz,Zygote.dropgrad(Δx),Zygote.dropgrad(Δy),Zygote.dropgrad(Δz),Zygote.dropgrad(Nx),Zygote.dropgrad(Ny),Zygote.dropgrad(Nz))
	# mag,mn = calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	mag,mn = calc_kpg(kz,g.g⃗)
	# ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	H,ω² = solve_ω²(kz,ε⁻¹,Δx,Δy,Δz;neigs,eigind,maxiter,tol)
	# println("ω² = $ω²")
	@show ω = sqrt(ω²)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	# ng = -ω / real( dot(Ha, kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,ds.mn), (2:4) ), ε⁻¹), (2:4)),ds.mn,ds.mag) ) )
	# ng = ω / H_Mₖ_H(Ha,ε⁻¹,mag,mn)
	ng = ω / real( dot(H, -vec( kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,mn), (2:4) ), ε⁻¹), (2:4)),mn,mag) ) ) )
	# ng = -ω / real( dot(Ha, kx_c2t( ifft( ε⁻¹_dot( fft( zx_t2c(Ha,Zygote.@showgrad(mn)), (2:4) ), ε⁻¹), (2:4)), Zygote.@showgrad(mn),Zygote.@showgrad(mag)) ) )
	( kz/ω, ng )
end

function solve_nω(kz::Array{<:Real},shapes::Vector{<:Shape},Δx,Δy,Δz,Nx,Ny,Nz;neigs=1,eigind=1,maxiter=3000,tol=1e-8)
	g::MaxwellGrid = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)  	# MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
	ε⁻¹::Array{Float64,5} = make_εₛ⁻¹(shapes,g)
	Hω = [solve_ω(kz[i],ε⁻¹,Δx,Δy,Δz;neigs,eigind,maxiter,tol) for i=1:length(kz)]
	ω² = [res[2] for res in Hω]
	ω = sqrt.(ω²)
	H = [res[1] for res in Hω]
	( kz ./ ω, [ ω[i] / H_Mₖ_H(H[i],ε⁻¹,calc_kpg(kz[i],Δx,Δy,Δz,Nx,Ny,Nz)...) for i=1:length(kz) ] ) # = (n, ng)
end

# using Zygote: @showgrad, dropgrad

# MkHa = Mₖ(Ha,ε⁻¹,kpg_mn,kpg_mag) #,ds.𝓕,ds.𝓕⁻¹)
# kxinds = [2; 1]
# kxscales = [-1.; 1.]
# @show size(H)
# temp = abs2.(H) #ε⁻¹_dot(zx_t2c(Ha,kpg_mn),ε⁻¹)
# Hastar = conj.(Ha)
# @tullio HMkH := Hastar[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * temp[a,i,j,k] * kpg_mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) nograd=(kxscales,kxinds,Hastar) fastmath=false verbose=2
# ng = ω / abs(HMkH)
# ng = sum(abs2,temp)
# ng = ω / real(H_Mₖ_H(H,ε⁻¹,kpg_mag,kpg_mn))
# ng = ω / H_Mₖ_H(H,ε⁻¹,Zygote.dropgrad(kpg_mag),Zygote.dropgrad(kpg_mn))

ext: jl


text:
@testset "SimpleGranularity" begin
    @test_throws ErrorException("Unknown type of simple granularity") SimpleGranularity("a")
end

@testset "DurationGranularity" begin
    @test_throws TypeError(:typeassert, "", String, 1) DurationGranularity(2, origin=1)
end

@testset "PeriodGranularity" begin
    @test_throws TypeError(:typeassert, "", String, 1) PeriodGranularity("a", origin=1)
    @test_throws TypeError(:typeassert, "", String, 1) PeriodGranularity("a", timezone=1)
end


text:
"""
	create_data_file(fname; data_dir, groups, datasets)

Create a data file with filename `fname` in directory `data_dir` and return the path.

Currently this always creates HDF5 files, for which optional inputs `groups` and `datasets`
specify names of group and dataset objects to create in the new HDF5 file.
"""
function create_data_file(fname; data_dir=joinpath(homedir(),"data"), groups=nothing,
	datasets=nothing)
	fpath 	= joinpath(data_dir, fname)
	fid 	= h5open(fpath, "cw")			# "cw"== "read-write, create if not existing"
	if !isnothing(groups)
		group_refs = [create_group(fid, gr_name) for gr_name in groups]
		ds_refs = [create_dataset(fid, ds_name) for ds_name in datasets]
	end
	close(fid)
	return fpath
end

function create_data_file(fname,datasets,sizes,types; data_dir=joinpath(homedir(),"data"), groups=nothing)
	fpath 	= joinpath(data_dir, fname)
	fid 	= h5open(fpath, "cw")			# "cw"== "read-write, create if not existing"
	if !isnothing(groups)
		group_refs = [create_group(fid, gr_name) for gr_name in groups]
        ds_refs = [create_dataset(fid, gr*"/"*ds, randn(tt,sz...)) for gr in groups, (ds,sz,tt) in zip(datasets,sizes,types)]
    else
        ds_refs = [create_dataset(fid, ds, randn(tt,sz...)) for (ds,sz,tt) in zip(datasets,sizes,types)]
    end
	close(fid)
	return fpath
end

function sweep_params(p_upper,p_lower,nps,ωs; fname=nothing, data_dir=joinpath(homedir(),"data"), grp_name=nothing)
	ranges = [range(p_lower[i],p_upper[i],length=nps[i]) for i=1:length(nps)]
	ps = collect(Iterators.product(ranges...)) # collect(Iterators.product(reverse(ranges)...))

	nω = length(ωs)
	fpath 	= joinpath(data_dir, fname)
	h5open(fpath, "cw") do fid			# "cw"== "read-write, create if not existing"
		ds_names = ["ngs_FH", "ngs_SH", "sumΔng²"]
		ds_sizes = [(nω,size(ps)...),(nω,size(ps)...),size(ps) ]
		ds_refs = [create_dataset(fid, nm, rand(sz...)) for (nm,sz) in zip(ds_names,ds_sizes)]
	end

	for pidx in CartesianIndices(ps)
		pp = ps[pidx]
		try
			sum_Δng²_FHSH(ωs,pp; pidx, fpath, grp_name)
		catch
			println("")
			println("pidx: $pidx")
			println("p: $pp")
			println("\tngs_FH: Error")
			println("\tngs_SH: Error")
			println("\tsum(Δng²): Error")
			println("")
		end
	end
end


text:
"""
    neighb_fun(price::Vector, demand, clear_error, pref, budget, capacity, ind_max::Int=40)

Return K neighbouring price vectors of the current vector p. The set is
the union of step search along the gradient of errors z and individual
course price adjustments.

# Parameters
- ind_max = the maxmimum number of individual price adjustments
"""
function neighb_fun(price::Vector, demand, clear_error, pref, budget, capacity, ind_max::Int=40)

    #gradient neighbours
    grad_neighb = grad_descent(price, clear_error)

    #individually adjusted neighbour prices
    if countnz(clear_error) <= ind_max
        indiv_neighb = [indiv_adjust(k, price, demand, clear_error, pref, budget, capacity) for k in 1:length(clear_error) if clear_error[k] != 0]
    else
        index  = sample(find(clear_error), ind_max, replace=false) # pick 40 random courses to adjust
        indiv_neighb = [indiv_adjust(k, price, demand, clear_error, pref, budget, capacity) for k in index]
    end

    # now get demand for all
    # neighbour_vect = [grad_neighb ; indiv_neighb]
    # D =d(neighbour_vect)
    # err =clearning_error(D[:course_demand])
    # return (neighbour_vect, D, err)

    return neighbour_vect = [grad_neighb ; indiv_neighb]
end

####### Gradient neighbours  #####

"Return the set of neighbours of price vector p found at different magnitudes
along the gradient of error. Increase the price of the most oversubscribed
course by magnitude, other course prices adjust relatively ."
function grad_descent(price::Vector, clear_error::Vector, steps=[10, 5, 1, 0.5, 0.1])
    relat_error = clear_error/maximum(clear_error) #relative clearing error
    grad_neighb = [max(price + (step * relat_error), 0) for step in steps] #price adjusted proportionaly to over/undersub, truncated at zero
    return grad_neighb
end

########## Individual adjustment ##########

"""
Adjust the price of course m to reduce demand by 1 if the course is oversubscribed
or drop the price to zero if the course is undersubscribed.
"""
function indiv_adjust(k::Int, price::Vector, demand, clear_error, pref, budget, capacity)
    adj_price = copy(price)

    if clear_error[k] < 0
        adj_price[k] = 0 # if course is undersubscribed, p = 0
    elseif clear_error[k] > 0
        adj_price[k] = decrease_demand(k, price, demand, clear_error, pref, budget, capacity)
    end

    return adj_price
end

"""
Return the adjusted price of course m such that its demand decreases by
exactly one.
"""
function decrease_demand(k::Int, price::Vector, demand, clear_error, pref, budget, capacity)
    S = length(pref) #number of students
    N = length(price) #number of courses
    min_increase = []

    for i in 1:S
        if dem["ind_demands"][i][k] == 1 #For all agents i demanding course k
        # 1) Find agent i's most prefered bundle not including x_m
            m1 = Model(solver=GurobiSolver())
            @variable(m1, x[1:N], Bin)
            @objective(m1, Max, x'*pref[i]*x) #max preferences
            @constraint(m1, sum(x) == capacity[i]) # course constraint
            @constraint(m1,  dot(price, x) <= budget[i]) # budget constraint
            @constraint(m1, x[k] == 0) # in set of bundles without course k
            status = solve(m1)
            o = getobjectivevalue(m1) # utility of most prefered bundle not including k

        # 2) Solve for minimum cost bundle in set of feasible solution with objective value greater than o
            m2 = Model(solver=GurobiSolver())
            @variable(m2, x[1:N], Bin)
            @objective(m2, Min, dot(price,x))
            @constraint(m2, x'*pref[i]*x >= o ) #bundles prefered to o
            @constraint(m2, sum(x) == capacity[i])
            @constraint(m2, x[k] == 1) #sechdule inclue course k
            status_2 = solve(m2)
            pi_opt = getobjectivevalue(m2) #cost of cheapest bundle including k

        # 3)Price increase for indiv i to drop course k
            if isnan(pi_opt) == false #drop infeasible results
                push!(min_increase, budget[i] - pi_opt + 1) #eps to be scaled to prices
            end
        end
    end
    return price[k] + minimum(min_increase)
end

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
module UnitfulOfu

using Unitful
import Unitful: N, m, W, g, Hz, @unit

#using Unitful: @unit
export @ofu_str

# OFU oz
# Now defined in Unitful
#@unit oz_ofu   "oz"        OFUOUNCE      (2*2834952//(10^5))*g     false

# OFU lbf
@unit lbf   "lbf"       OFUPound      (4.45)*N                  true

# OFU hp
@unit hp    "hp"        OFUHP         (735.4990027920646)*W     false

# OFU rpm
@unit rpm   "rpm"       OFURPM        (1//60)*Hz                false

# Define Nm
@unit Nm    "Nm"        OFUNm         1m*1N                     true

include("ofumacro.jl")

# Some gymnastics required here because if we precompile, we cannot add to
# Unitful.basefactors at compile time and expect the changes to persist to runtime.
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulOfu)
end

end # module


text:
module _Algorithm

using Brainstorm._Algorithm
using Base.Test

include("binarysearch.jl")

function test_all()
    println("")

    test_binarysearch_all()
end

end


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

## type join (closest common ancestor, or least upper bound) ##

"""
    typejoin(T, S)

Return the closest common ancestor of `T` and `S`, i.e. the narrowest type from which
they both inherit.
"""
typejoin() = (@_pure_meta; Bottom)
typejoin(@nospecialize(t)) = (@_pure_meta; t)
typejoin(@nospecialize(t), ts...) = (@_pure_meta; typejoin(t, typejoin(ts...)))
typejoin(@nospecialize(a), @nospecialize(b)) = (@_pure_meta; join_types(a, b, typejoin))

function join_types(@nospecialize(a), @nospecialize(b), f::Function)
    if a <: b
        return b
    elseif b <: a
        return a
    elseif isa(a,UnionAll)
        return UnionAll(a.var, join_types(a.body, b, f))
    elseif isa(b,UnionAll)
        return UnionAll(b.var, join_types(a, b.body, f))
    elseif isa(a,TypeVar)
        return f(a.ub, b)
    elseif isa(b,TypeVar)
        return f(a, b.ub)
    elseif isa(a,Union)
        a′ = f(a.a,a.b)
        return a′ === a ? typejoin(a, b) : f(a′, b)
    elseif isa(b,Union)
        b′ = f(b.a,b.b)
        return b′ === b ? typejoin(a, b) : f(a, b′)
    elseif a <: Tuple
        if !(b <: Tuple)
            return Any
        end
        ap, bp = a.parameters, b.parameters
        lar = length(ap)::Int; lbr = length(bp)::Int
        if lar == 0
            return Tuple{Vararg{tailjoin(bp,1,f)}}
        end
        if lbr == 0
            return Tuple{Vararg{tailjoin(ap,1,f)}}
        end
        laf, afixed = full_va_len(ap)
        lbf, bfixed = full_va_len(bp)
        if laf < lbf
            if isvarargtype(ap[lar]) && !afixed
                c = Vector{Any}(uninitialized, laf)
                c[laf] = Vararg{f(unwrapva(ap[lar]), tailjoin(bp,laf,f))}
                n = laf-1
            else
                c = Vector{Any}(uninitialized, laf+1)
                c[laf+1] = Vararg{tailjoin(bp,laf+1,f)}
                n = laf
            end
        elseif lbf < laf
            if isvarargtype(bp[lbr]) && !bfixed
                c = Vector{Any}(uninitialized, lbf)
                c[lbf] = Vararg{f(unwrapva(bp[lbr]), tailjoin(ap,lbf,f))}
                n = lbf-1
            else
                c = Vector{Any}(uninitialized, lbf+1)
                c[lbf+1] = Vararg{tailjoin(ap,lbf+1,f)}
                n = lbf
            end
        else
            c = Vector{Any}(uninitialized, laf)
            n = laf
        end
        for i = 1:n
            ai = ap[min(i,lar)]; bi = bp[min(i,lbr)]
            ci = f(unwrapva(ai),unwrapva(bi))
            c[i] = i == length(c) && (isvarargtype(ai) || isvarargtype(bi)) ? Vararg{ci} : ci
        end
        return Tuple{c...}
    elseif b <: Tuple
        return Any
    end
    while b !== Any
        if a <: b.name.wrapper
            while a.name !== b.name
                a = supertype(a)
            end
            aprimary = unwrap_unionall(a.name.wrapper)
            # join on parameters
            n = length(a.parameters)
            if n == 0
                return aprimary
            end
            p = Vector{Any}(uninitialized, n)
            for i = 1:n
                ai, bi = a.parameters[i], b.parameters[i]
                if ai === bi || (isa(ai,Type) && isa(bi,Type) && typeseq(ai,bi))
                    p[i] = ai
                else
                    p[i] = aprimary.parameters[i]
                end
            end
            return rewrap_unionall(a.name.wrapper{p...}, a.name.wrapper)
        end
        b = supertype(b)
    end
    return Any
end

"""
    promote_typejoin(T, S)

Compute a type that contains both `T` and `S`, which could be
either a parent of both types, or a `Union` if appropriate.
Falls back to [`typejoin`](@ref).
"""
promote_typejoin(@nospecialize(a), @nospecialize(b)) =
    (@_pure_meta; join_types(a, b, promote_typejoin))
promote_typejoin(::Type{Nothing}, ::Type{T}) where {T} =
    isconcretetype(T) ? Union{T, Nothing} : Any
promote_typejoin(::Type{T}, ::Type{Nothing}) where {T} =
    isconcretetype(T) ? Union{T, Nothing} : Any
promote_typejoin(::Type{Missing}, ::Type{T}) where {T} =
    isconcretetype(T) ? Union{T, Missing} : Any
promote_typejoin(::Type{T}, ::Type{Missing}) where {T} =
    isconcretetype(T) ? Union{T, Missing} : Any
promote_typejoin(::Type{Nothing}, ::Type{Missing}) = Union{Nothing, Missing}
promote_typejoin(::Type{Missing}, ::Type{Nothing}) = Union{Nothing, Missing}
promote_typejoin(::Type{Nothing}, ::Type{Nothing}) = Nothing
promote_typejoin(::Type{Missing}, ::Type{Missing}) = Missing

# Returns length, isfixed
function full_va_len(p)
    isempty(p) && return 0, true
    last = p[end]
    if isvarargtype(last)
        N = unwrap_unionall(last).parameters[2]
        if isa(N, Integer)
            return (length(p) + N - 1)::Int, true
        end
        return length(p)::Int, false
    end
    return length(p)::Int, true
end

# reduce join_types over A[i:end]
function tailjoin(A, i, f::Function)
    if i > length(A)
        return unwrapva(A[end])
    end
    t = Bottom
    for j = i:length(A)
        t = f(t, unwrapva(A[j]))
    end
    return t
end

## promotion mechanism ##

"""
    promote_type(type1, type2)

Promotion refers to converting values of mixed types to a single common type.
`promote_type` represents the default promotion behavior in Julia when
operators (usually mathematical) are given arguments of differing types.
`promote_type` generally tries to return a type which can at least approximate
most values of either input type without excessively widening.  Some loss is
tolerated; for example, `promote_type(Int64, Float64)` returns
[`Float64`](@ref) even though strictly, not all [`Int64`](@ref) values can be
represented exactly as `Float64` values.

```jldoctest
julia> promote_type(Int64, Float64)
Float64

julia> promote_type(Int32, Int64)
Int64

julia> promote_type(Float32, BigInt)
BigFloat

julia> promote_type(Int16, Float16)
Float16

julia> promote_type(Int64, Float16)
Float16

julia> promote_type(Int8, UInt16)
UInt16
```
"""
function promote_type end

promote_type()  = Bottom
promote_type(T) = T
promote_type(T, S, U, V...) = (@_inline_meta; promote_type(T, promote_type(S, U, V...)))

promote_type(::Type{Bottom}, ::Type{Bottom}) = Bottom
promote_type(::Type{T}, ::Type{T}) where {T} = T
promote_type(::Type{T}, ::Type{Bottom}) where {T} = T
promote_type(::Type{Bottom}, ::Type{T}) where {T} = T

function promote_type(::Type{T}, ::Type{S}) where {T,S}
    @_inline_meta
    # Try promote_rule in both orders. Typically only one is defined,
    # and there is a fallback returning Bottom below, so the common case is
    #   promote_type(T, S) =>
    #   promote_result(T, S, result, Bottom) =>
    #   typejoin(result, Bottom) => result
    promote_result(T, S, promote_rule(T,S), promote_rule(S,T))
end

"""
    promote_rule(type1, type2)

Specifies what type should be used by [`promote`](@ref) when given values of types `type1` and
`type2`. This function should not be called directly, but should have definitions added to
it for new types as appropriate.
"""
function promote_rule end

promote_rule(::Type{<:Any}, ::Type{<:Any}) = Bottom
# To fix ambiguities
promote_rule(::Type{Any}, ::Type{<:Any}) = Any
promote_rule(::Type{<:Any}, ::Type{Any}) = Any
promote_rule(::Type{Any}, ::Type{Any}) = Any

promote_result(::Type{<:Any},::Type{<:Any},::Type{T},::Type{S}) where {T,S} = (@_inline_meta; promote_type(T,S))
# If no promote_rule is defined, both directions give Bottom. In that
# case use typejoin on the original types instead.
promote_result(::Type{T},::Type{S},::Type{Bottom},::Type{Bottom}) where {T,S} = (@_inline_meta; typejoin(T, S))

"""
    promote(xs...)

Convert all arguments to a common type, and return them all (as a tuple).
If no arguments can be converted, an error is raised.

# Examples
```jldoctest
julia> promote(Int8(1), Float16(4.5), Float32(4.1))
(1.0f0, 4.5f0, 4.1f0)
```
"""
function promote end

function _promote(x::T, y::S) where {T,S}
    @_inline_meta
    R = promote_type(T, S)
    return (convert(R, x), convert(R, y))
end
promote_typeof(x) = typeof(x)
promote_typeof(x, xs...) = (@_inline_meta; promote_type(typeof(x), promote_typeof(xs...)))
function _promote(x, y, z)
    @_inline_meta
    R = promote_typeof(x, y, z)
    return (convert(R, x), convert(R, y), convert(R, z))
end
function _promote(x, y, zs...)
    @_inline_meta
    R = promote_typeof(x, y, zs...)
    return (convert(R, x), convert(R, y), convert(Tuple{Vararg{R}}, zs)...)
end
# TODO: promote(x::T, ys::T...) where {T} here to catch all circularities?

## promotions in arithmetic, etc. ##

# Because of the promoting fallback definitions for Number, we need
# a special case for undefined promote_rule on numeric types.
# Otherwise, typejoin(T,S) is called (returning Number) so no conversion
# happens, and +(promote(x,y)...) is called again, causing a stack
# overflow.
function promote_result(::Type{T},::Type{S},::Type{Bottom},::Type{Bottom}) where {T<:Number,S<:Number}
    @_inline_meta
    promote_to_supertype(T, S, typejoin(T,S))
end

# promote numeric types T and S to typejoin(T,S) if T<:S or S<:T
# for example this makes promote_type(Integer,Real) == Real without
# promoting arbitrary pairs of numeric types to Number.
promote_to_supertype(::Type{T}, ::Type{T}, ::Type{T}) where {T<:Number}           = (@_inline_meta; T)
promote_to_supertype(::Type{T}, ::Type{S}, ::Type{T}) where {T<:Number,S<:Number} = (@_inline_meta; T)
promote_to_supertype(::Type{T}, ::Type{S}, ::Type{S}) where {T<:Number,S<:Number} = (@_inline_meta; S)
promote_to_supertype(::Type{T}, ::Type{S}, ::Type) where {T<:Number,S<:Number} =
    error("no promotion exists for ", T, " and ", S)

promote() = ()
promote(x) = (x,)

function promote(x, y)
    @_inline_meta
    px, py = _promote(x, y)
    not_sametype((x,y), (px,py))
    px, py
end
function promote(x, y, z)
    @_inline_meta
    px, py, pz = _promote(x, y, z)
    not_sametype((x,y,z), (px,py,pz))
    px, py, pz
end
function promote(x, y, z, a...)
    p = _promote(x, y, z, a...)
    not_sametype((x, y, z, a...), p)
    p
end

promote(x::T, y::T, zs::T...) where {T} = (x, y, zs...)

not_sametype(x::T, y::T) where {T} = sametype_error(x)

not_sametype(x, y) = nothing

function sametype_error(input)
    @_noinline_meta
    error("promotion of types ",
          join(map(x->string(typeof(x)), input), ", ", " and "),
          " failed to change any arguments")
end

+(x::Number, y::Number) = +(promote(x,y)...)
*(x::Number, y::Number) = *(promote(x,y)...)
-(x::Number, y::Number) = -(promote(x,y)...)
/(x::Number, y::Number) = /(promote(x,y)...)

"""
    ^(x, y)

Exponentiation operator. If `x` is a matrix, computes matrix exponentiation.

If `y` is an `Int` literal (e.g. `2` in `x^2` or `-3` in `x^-3`), the Julia code
`x^y` is transformed by the compiler to `Base.literal_pow(^, x, Val(y))`, to
enable compile-time specialization on the value of the exponent.
(As a default fallback we have `Base.literal_pow(^, x, Val(y)) = ^(x,y)`,
where usually `^ == Base.^` unless `^` has been defined in the calling
namespace.)

```jldoctest
julia> 3^5
243

julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> A^3
2×2 Array{Int64,2}:
 37   54
 81  118
```
"""
^(x::Number, y::Number) = ^(promote(x,y)...)

fma(x::Number, y::Number, z::Number) = fma(promote(x,y,z)...)
muladd(x::Number, y::Number, z::Number) = muladd(promote(x,y,z)...)

==(x::Number, y::Number) = (==)(promote(x,y)...)
<( x::Real, y::Real)     = (< )(promote(x,y)...)
<=(x::Real, y::Real)     = (<=)(promote(x,y)...)

div(x::Real, y::Real) = div(promote(x,y)...)
fld(x::Real, y::Real) = fld(promote(x,y)...)
cld(x::Real, y::Real) = cld(promote(x,y)...)
rem(x::Real, y::Real) = rem(promote(x,y)...)
mod(x::Real, y::Real) = mod(promote(x,y)...)

mod1(x::Real, y::Real) = mod1(promote(x,y)...)
fld1(x::Real, y::Real) = fld1(promote(x,y)...)

max(x::Real, y::Real) = max(promote(x,y)...)
min(x::Real, y::Real) = min(promote(x,y)...)
minmax(x::Real, y::Real) = minmax(promote(x, y)...)

# "Promotion" that takes a function into account and tries to preserve
# non-concrete types. These are meant to be used mainly by elementwise
# operations, so it is advised against overriding them
_default_type(T::Type) = (@_inline_meta; T)

if isdefined(Core, :Compiler)
    const _return_type = Core.Compiler.return_type
else
    _return_type(@nospecialize(f), @nospecialize(t)) = Any
end

"""
    promote_op(f, argtypes...)

Guess what an appropriate container eltype would be for storing results of
`f(::argtypes...)`. The guess is in part based on type inference, so can change any time.

!!! warning
    In pathological cases, the type returned by `promote_op(f, argtypes...)` may not even
    be a supertype of the return value of `f(::argtypes...)`. Therefore, `promote_op`
    should _not_ be used e.g. in the preallocation of an output array.

!!! warning
    Due to its fragility, use of `promote_op` should be avoided. It is preferable to base
    the container eltype on the type of the actual elements. Only in the absence of any
    elements (for an empty result container), it may be unavoidable to call `promote_op`.
"""
promote_op(::Any...) = (@_inline_meta; Any)
function promote_op(f, ::Type{S}) where S
    @_inline_meta
    TT = Tuple{_default_type(S)}
    T = _return_type(f, TT)
    isdispatchtuple(Tuple{S}) && return isdispatchtuple(Tuple{T}) ? T : Any
    return typejoin(S, T)
end
function promote_op(f, ::Type{R}, ::Type{S}) where {R,S}
    @_inline_meta
    TT = Tuple{_default_type(R), _default_type(S)}
    T = _return_type(f, TT)
    isdispatchtuple(Tuple{R}) && isdispatchtuple(Tuple{S}) && return isdispatchtuple(Tuple{T}) ? T : Any
    return typejoin(R, S, T)
end

## catch-alls to prevent infinite recursion when definitions are missing ##

no_op_err(name, T) = error(name," not defined for ",T)
(+)(x::T, y::T) where {T<:Number} = no_op_err("+", T)
(*)(x::T, y::T) where {T<:Number} = no_op_err("*", T)
(-)(x::T, y::T) where {T<:Number} = no_op_err("-", T)
(/)(x::T, y::T) where {T<:Number} = no_op_err("/", T)
(^)(x::T, y::T) where {T<:Number} = no_op_err("^", T)

fma(x::T, y::T, z::T) where {T<:Number} = no_op_err("fma", T)
fma(x::Integer, y::Integer, z::Integer) = x*y+z
muladd(x::T, y::T, z::T) where {T<:Number} = x*y+z

(&)(x::T, y::T) where {T<:Integer} = no_op_err("&", T)
(|)(x::T, y::T) where {T<:Integer} = no_op_err("|", T)
xor(x::T, y::T) where {T<:Integer} = no_op_err("xor", T)

(==)(x::T, y::T) where {T<:Number} = x === y
(< )(x::T, y::T) where {T<:Real} = no_op_err("<" , T)
(<=)(x::T, y::T) where {T<:Real} = no_op_err("<=", T)

rem(x::T, y::T) where {T<:Real} = no_op_err("rem", T)
mod(x::T, y::T) where {T<:Real} = no_op_err("mod", T)

min(x::Real) = x
max(x::Real) = x
minmax(x::Real) = (x, x)

max(x::T, y::T) where {T<:Real} = select_value(y < x, x, y)
min(x::T, y::T) where {T<:Real} = select_value(y < x, y, x)
minmax(x::T, y::T) where {T<:Real} = y < x ? (y, x) : (x, y)

flipsign(x::T, y::T) where {T<:Signed} = no_op_err("flipsign", T)

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
abstract type AbstractFVAlgorithm end

function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::AbstractFVAlgorithm, noscalar::Bool, ::Type{Val{false}})
    # This works assuming uniform mesh TODO: trait for uniform meshes
    dx = cell_volume(mesh, 1)
    #update vector
    for j in node_indices(mesh)
        if noscalar
            fluxes[:,j] .= update_flux_value(u,j,dt,dx,mesh,Flux,alg)
        else
            fluxes[j] = update_flux_value(u,j,dt,dx,mesh,Flux,alg)
        end
    end
end

function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::AbstractFVAlgorithm, noscalar::Bool, ::Type{Val{true}})
    # This works assuming uniform mesh TODO: trait for uniform meshes
    dx = cell_volume(mesh, 1)
    #update vector
    Threads.@threads for j in node_indices(mesh)
        if noscalar
            fluxes[:,j] .= update_flux_value(u,j,dt,dx,mesh,Flux,alg)
        else
            fluxes[j] = update_flux_value(u,j,dt,dx,mesh,Flux,alg)
        end
    end
end


text:
"""
    getthreads(token::String)
Takes in the Slack App Auth token, channel ID, and Time Signature of the parent message.
Return the first 50 replies ina  thread.

https://api.slack.com/methods/conversations.replies#arg_ts
"""
function getthreads(token::String, channelid::String, ts::String)
    HTTP.setuseragent!("HTTP.jl")

    # TODO: 50 is an arbitary number, see if there is a limit.
    url = "https://slack.com/api/conversations.replies?token=$(token)&channel=$(channelid)&ts=$(ts)&limit=50"

    r = HTTP.get(url, ["Content-Type" => "application/json"])
    message_holder = []

    json_obj = JSON.parse(String(r.body))
    message_holder = []
    for (k,v) in json_obj
        if occursin("messages", k)
            for item in v

                # Declare the holder variables
                edited = false
                ts = ""
                thread_ts = ""
                user = ""
                text = ""

                # println("Starting to look through next Dict/Object") #DEBUG
                for (key, value) in item
                    if key == "user"
                        user = value

                    elseif key == "edited"
                        edited = true

                    elseif key == "ts"
                        ts = value

                    elseif key == "thread_ts"
                        thread_ts = value

                    elseif key == "text"
                        text = value
                    end
                end

                # Create the message object.. Handle a failed parse and object creation error.
                try
                    new_tm = threadedMessage(edited, ts, thread_ts, user, text)

                    # Push it into the holder array
                    push!(message_holder, new_tm)
                catch
                    println("Enabled to created a threaded message with the following details: $(edited), $(ts), $(thread_ts), $(user), $(text)")
                end

            end
        end
    end

    # Create and Return the Thread object (just an array of messages)
    return(message_holder)

end


text:
# This file was generated, do not modify it. # hide
caravan  = dataset("ISLR", "Caravan")
size(caravan)


text:
export train_pq, quantize_pq

"""
    quantize_pq(X, C, V=false) -> B

Given data and PQ codeboks, quantize.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `C::Vector{Matrix{T}}`: `m`-long vector with `d/m`-by-`h` matrix entries. Each matrix is a PQ codebook.
- `V::Bool`: Whether to print progress

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes that approximate `X`
"""
function quantize_pq(
  X::Matrix{T},             # d-by-n matrix of data points to quantize
  C::Vector{Matrix{T}},     # m-long array. Each entry is a d-by-h codebooks
  V::Bool=false) where T <: AbstractFloat

  d, n = size(X)
  m    = length(C)
  h    = size(C[1], 2)
  B    = Vector{Vector{Int}}(undef, m) # codes
  for i = 1:m; B[i] = zeros(Int, n) end # Allocate codes

  subdims = splitarray(1:d, m)

  # auxiliary variables for update_assignments! function
  costs     = zeros(Float32, n)
  counts    = zeros(Int, h)
  to_update = zeros(Bool, h)
  unused    = Int[]

  for i = 1:m
    if V print("Encoding on codebook $i / $m... ") end
    # Find distances from X to codebook
    dmat = Distances.pairwise(Distances.SqEuclidean(), C[i], X[subdims[i],:])
    Clustering.update_assignments!(dmat, true, B[i], costs, counts, to_update, unused)
    if V println("done"); end
  end

  B = hcat(B...)
  B = convert(Matrix{Int16}, B)
  collect(B')
end

"""
    train_pq(X, m, h, niter=25, V=false) -> C, B, error

Trains a product quantizer.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `m::Integer`: Number of codebooks
- `h::Integer`: Number of entries in each codebook (typically 256)
- `niter::Integer`: Number of k-means iterations to use
- `V::Bool`: Whether to print progress

# Returns
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook of size approximately `d/m`-by-`h`.
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `error::T`: The quantization error after training
"""
function train_pq(
  X::Matrix{T},         # d-by-n matrix of data points to train on.
  m::Integer,           # number of codebooks
  h::Integer,           # number of entries per codebook
  niter::Integer=25,    # number of optimization iterations
  V::Bool=false) where T <: AbstractFloat

  d, n = size(X)

  C = Vector{Matrix{T}}(undef, m) # codebooks

  B       = zeros(Int16, n, m) # codes
  subdims = splitarray(1:d, m) # subspaces

  for i = 1:m
    if V print("Working on codebook $i / $m... "); end
    # FAISS uses 25 iterations by default
    # https://github.com/facebookresearch/faiss/blob/master/Clustering.cpp#L28
    cluster = kmeans(X[ subdims[i],: ], h, init=:kmpp, maxiter=niter)
    C[i], B[:,i] = cluster.centers, cluster.assignments

    if V
      println("done.")
      println("  Ran for $(cluster.iterations) iterations")
      println("  Error in subspace is $(cluster.totalcost ./ n)")
      println("  Converged: $(cluster.converged)")
    end
  end
  B = collect(B')
  error = qerror_pq(X, B, C)
  return C, B, error
end

"""
  Run a complete experiment: train, encode base and compute recall.
"""
function experiment_pq(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === Train ===
  C, B, train_error = train_pq(Xt, m, h, niter, V)
  if V; @printf("Error in training is %e\n", train_error); end

  # === Encode the base set ===
  B_base     = quantize_pq(Xb, C, V)
  base_error = qerror_pq(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # === Compute recall ===
  if V; println("Querying m=$m ... "); end
  b = Int(log2(h) * m)
  @time dists, idx = linscan_pq(B_base, Xq, C, b, knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, B_base, recall
end

"""
  Run a complete experiment where the training set and base set are the same.
"""
function experiment_pq_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === Train ===
  C, B, train_error = train_pq(Xt, m, h, niter, V)
  if V; @printf("Error in training is %e\n", train_error); end

  # === Compute recall ===
  if V; println("Querying m=$m ... "); end
  b = Int(log2(h) * m)
  @time dists, idx = linscan_pq(B, Xq, C, b, knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, recall
end


text:
#HadamardProd

export HadamardProd

"""
`HadamardProd(A::AbstractOperator,B::AbstractOperator)`

Create an operator `P` such that:

`P*x == (Ax).*(Bx)`

# Example

```julia
julia> A,B = Sin(3), Cos(3);

julia> P = HadamardProd(A,B)
sin.*cos  ℝ^3 -> ℝ^3

julia> x = randn(3);

julia> P*x == (sin.(x).*cos.(x))
true

```
"""
struct HadamardProd{
                 L1 <: AbstractOperator,
                 L2 <: AbstractOperator,
                 C <: AbstractArray,
                 D <: AbstractArray,
                } <: NonLinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufD::D
  function HadamardProd(A::L1, B::L2, bufA::C, bufB::C, bufD::D) where {L1,L2,C,D}
    if size(A) != size(B)
      throw(DimensionMismatch("Cannot compose operators"))
    end
    new{L1,L2,C,D}(A,B,bufA,bufB,bufD)
  end
end

struct HadamardProdJac{
                    L1 <: AbstractOperator,
                    L2 <: AbstractOperator,
                    C <: AbstractArray,
                    D <: AbstractArray,
                   } <: LinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufD::D
end

# Constructors
function HadamardProd(A::AbstractOperator,B::AbstractOperator)
  s,t = size(A,1), codomainType(A)
  bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(B,1), codomainType(B)
  bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(A,2), domainType(A)
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  HadamardProd(A,B,bufA,bufB,bufD)
end

# Jacobian
function Jacobian(P::HadamardProd{L1,L2,C,D}, x::AbstractArray) where {L1,L2,C,D}
  JA, JB = Jacobian(P.A, x), Jacobian(P.B, x)
  HadamardProdJac{typeof(JA),typeof(JB),C,D}(JA,JB,P.bufA,P.bufB,P.bufD)
end

# Mappings
function mul!(y, P::HadamardProd{L1,L2,C,D}, b) where {L1,L2,C,D}
  mul!(P.bufA,P.A,b)
  mul!(P.bufB,P.B,b)
  y .= P.bufA .* P.bufB
  return y
end

function mul!(y, J::AdjointOperator{HadamardProdJac{L1,L2,C,D}}, b) where {L1,L2,C,D}
  #y .= J.A.B' * ( J.A.bufA .*b ) + J.A.A' * ( J.A.bufB .* b )
  J.A.bufA .*= b
  mul!(y, J.A.B', J.A.bufA)
  J.A.bufB .*= b
  mul!(J.A.bufD, J.A.A', J.A.bufB)
  y .+= J.A.bufD
  return y
end

size(P::Union{HadamardProd,HadamardProdJac}) = (size(P.A,1),size(P.A,2))

fun_name(L::Union{HadamardProd,HadamardProdJac}) = fun_name(L.A)*".*"*fun_name(L.B)

domainType(L::Union{HadamardProd,HadamardProdJac})   = domainType(L.A)
codomainType(L::Union{HadamardProd,HadamardProdJac}) = codomainType(L.A)

# utils
function permute(P::HadamardProd{L1,L2,C,D},
                 p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition}
  HadamardProd(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,ArrayPartition(P.bufD.x[p]) )
end

remove_displacement(P::HadamardProd) =
HadamardProd(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufD)


text:
using QuantumStatePlots
using QuantumStateBase
using Test

@testset "QuantumStatePlots.jl" begin

ENV["GKSwstype"]="nul"

@testset "plot wigner" begin
    x_range = -5:1.0:5
    p_range = -5:1.0:5
    wf = WignerFunction(x_range, p_range)

    state = VacuumState()
    ws = wf(state)

    file_path = "wigner.png"

    plot_wigner(ws, Heatmap, file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)

    plot_wigner(ws, Contour, file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)

    plot_wigner(ws, Surface, file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)
end

@testset "plot ρ" begin
    file_path = "rho.png"

    plot_ρ(VacuumState(), file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)

    plot_ρ(VacuumState(), state_n=5, file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)
end

@testset "plot all" begin
    x_range = -5:1.0:5
    p_range = -5:1.0:5
    wf = WignerFunction(x_range, p_range)

    state = VacuumState()
    ws = wf(state)

    file_path = "all.png"

    plot_all(ws, state, file_path=file_path)
    @test isfile(file_path)
    isfile(file_path) && rm(file_path)
end

end


text:
__precompile__()

module Graphics

using Compat
using Compat.LinearAlgebra
using Colors
using NaNMath

import Base: +, -, *, /, &, fill
import Compat.LinearAlgebra: norm

"""
Graphics defines an API for drawing in two dimensions.

- Geometric primitives: `Vec2`, `Point`, `BoundingBox`

- Geometry API: `aspect_ratio`, `center`, `deform`, `diagonal`,
  `isinside`, `shift`, `height`, `width`, `xmin`, `xmax`, `ymin`,
  `ymax`, `xrange`, `yrange`

- 2d drawing contexts: `GraphicsDevice`, `GraphicsContext`, `creategc`, `getgc`

- Coordinate systems: `set_coordinates`, `reset_transform`, `rotate`,
  `scale`, `translate`, `user_to_device!`, `device_to_user!`,
  `user_to_device_distance!`, `device_to_user_distance!`,
  `user_to_device`, `device_to_user`

- Lines: `set_line_width`, `set_dash`

- Colors and painting (drawing attributes): `set_source`,
  `set_source_rgb`, `set_source_rgba`, `save`, `restore`

- Clipping: `clip`, `clip_preserve`, `reset_clip`, `inner_canvas`

- Paths: `move_to`, `line_to`, `rel_line_to`, `rel_move_to`,
  `new_path`, `new_sub_path`, `close_path`, `arc`

- High-level paths: `rectangle`, `circle`, `polygon`

- Fill and stroke: `fill`, `fill_preserve`, `paint`, `stroke`,
  `stroke_preserve`, `stroke_transformed`,
  `stroke_transformed_preserve`
"""
Graphics

export
    # Part 1. 2D Geometry
    Vec2, Point, BoundingBox,
    # limits in world coordinates
    isinside, xmin, xmax, ymin, ymax, center, xrange, yrange,
    aspect_ratio, with_aspect_ratio, diagonal, shift, deform,
    # TODO: more

    # TODO: 3D geometry

    # Part 2. 2D Drawing
    # device and context
    GraphicsDevice, GraphicsContext, creategc, getgc,
    # width, height are in user (world) coordinates for geometric objects,
    # but in device coordinates for GraphicsDevice, GraphicsContext, and
    # other concrete things like windows and widgets.
    width, height,

    # drawing attribute manipulation
    save, restore, set_line_width, set_dash, set_source_rgb, set_source_rgba,
    set_source,

    # coordinate systems
    reset_transform, set_coordinates, rotate, scale, translate, user_to_device!,
    device_to_user!, user_to_device_distance!, device_to_user_distance!,
    user_to_device, device_to_user,

    # clipping
    clip, clip_preserve, reset_clip, inner_canvas,

    # path primitives
    move_to, line_to, rel_line_to, rel_move_to, new_path, new_sub_path,
    close_path, arc,

    # fill and stroke
    fill, fill_preserve, paint, stroke, stroke_preserve,
    stroke_transformed, stroke_transformed_preserve,

    # derived path operations
    rectangle, circle, polygon

    # TODO: text drawing API

    # TODO: rendering pipeline API

# Part 1. geometric primitives

"""
    Vec2(x, y) -> v

Create a Cartesian representation `v` of a vector (or point) in two dimensions.
"""
struct Vec2
    x::Float64
    y::Float64
end

"""
    Point(x, y) -> p

Create a Cartesian representation `p` of a point in two dimensions.
"""
const Point = Vec2

(+)(a::Vec2, b::Vec2) = Vec2(a.x + b.x, a.y + b.y)
(-)(a::Vec2, b::Vec2) = Vec2(a.x - b.x, a.y - b.y)
(*)(p::Vec2, s::Real) = Vec2(p.x*s, p.y*s)
(/)(p::Vec2, s::Real) = Vec2(p.x/s, p.y/s)
(*)(s::Real, p::Vec2) = p*s

"""
    rotate(p::Vec2, angle::Real, o::Vec2) -> pnew

Rotate `p` around `o` by `angle`.
"""
function rotate(p::Vec2, angle::Real, o::Vec2)
    c = cos(angle)
    s = sin(angle)
    d = p - o
    Vec2(o.x + c*d.x - s*d.y, o.y + s*d.x + c*d.y)
end
rotate(p::Vec2, angle::Real) = rotate(p, angle, Vec2(0.,0.))

norm(p::Vec2) = hypot(p.x, p.y)

"""
    BoundingBox(xmin, xmax, ymin, ymax) -> bb

Create a representation `bb` of a rectangular region, specifying the
coordinates of the horizontal (x) and vertical (y) edges.
"""
struct BoundingBox
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
end

BoundingBox() = BoundingBox(NaN, NaN, NaN, NaN)

"""
    BoundingBox(p0::Point, points::Point...) -> bb

Compute the BoundingBox `bb` that minimally encloses all of the input points.
"""
function BoundingBox(p0::Point, points::Point...)
    xmin, xmax, ymin, ymax = p0.x, p0.x, p0.y, p0.y
    for p in points
        xmin = NaNMath.min(xmin, p.x)
        xmax = NaNMath.max(xmax, p.x)
        ymin = NaNMath.min(ymin, p.y)
        ymax = NaNMath.max(ymax, p.y)
    end
    return BoundingBox(xmin, xmax, ymin, ymax)
end

"""
    BoundingBox(bb0::BoundingBox, bboxes::BoundingBox...) -> bb

Compute the BoundingBox `bb` that minimally encloses all of the input boxes.
"""
function BoundingBox(bb0::BoundingBox, bboxes::BoundingBox...)
    xmin, xmax, ymin, ymax = bb0.xmin, bb0.xmax, bb0.ymin, bb0.ymax
    for bb in bboxes
        xmin = NaNMath.min(xmin, bb.xmin)
        xmax = NaNMath.max(xmax, bb.xmax)
        ymin = NaNMath.min(ymin, bb.ymin)
        ymax = NaNMath.max(ymax, bb.ymax)
    end
    return BoundingBox(xmin, xmax, ymin, ymax)
end

width(bb::BoundingBox) = bb.xmax - bb.xmin
height(bb::BoundingBox) = bb.ymax - bb.ymin
diagonal(bb) = hypot(width(bb), height(bb))
"""
    aspect_ratio(bb::BoundingBox) -> r

Compute the ratio `r` of the height and width of `bb`.
"""
aspect_ratio(bb) = height(bb)/width(bb)

xmin(bb::BoundingBox) = bb.xmin
xmax(bb::BoundingBox) = bb.xmax
ymin(bb::BoundingBox) = bb.ymin
ymax(bb::BoundingBox) = bb.ymax

center(x) = Point((xmin(x)+xmax(x))/2, (ymin(x)+ymax(x))/2)
xrange(x) = xmin(x), xmax(x)
yrange(x) = ymin(x), ymax(x)

"""
    bb1 + bb2 -> bb

Compute the BoundingBox `bb` that minimally contains `bb1` and `bb2`
"""
function (+)(bb1::BoundingBox, bb2::BoundingBox)
    BoundingBox(NaNMath.min(bb1.xmin, bb2.xmin),
                NaNMath.max(bb1.xmax, bb2.xmax),
                NaNMath.min(bb1.ymin, bb2.ymin),
                NaNMath.max(bb1.ymax, bb2.ymax))
end

"""
    bb1 & bb2 -> bb

Compute the intersection of two BoundingBoxes.
"""
function (&)(bb1::BoundingBox, bb2::BoundingBox)
    BoundingBox(NaNMath.max(bb1.xmin, bb2.xmin),
                NaNMath.min(bb1.xmax, bb2.xmax),
                NaNMath.max(bb1.ymin, bb2.ymin),
                NaNMath.min(bb1.ymax, bb2.ymax))
end

"""
    deform(bb::BoundingBox, Δl, Δr, Δt, Δb) -> bbnew

Add `Δl` (left), `Δr` (right), `Δt` (top), and `Δb` (bottom) to the
edges of a BoundingBox.
"""
function deform(bb::BoundingBox, dl, dr, dt, db)
    BoundingBox(bb.xmin + dl, bb.xmax + dr, bb.ymin + dt, bb.ymax + db)
end

"""
    shift(bb::BoundingBox, Δx, Δy) -> bbnew

Shift center by (Δx,Δy), keeping width & height fixed.
"""
function shift(bb::BoundingBox, dx, dy)
    BoundingBox(bb.xmin + dx, bb.xmax + dx, bb.ymin + dy, bb.ymax + dy)
end

"""
    s*bb -> bbnew
    bb*s -> bbnew

Scale width & height of BoundingBox `bb` by `s`, keeping center fixed.
"""
function (*)(bb::BoundingBox, s::Real)
    dw = 0.5*(s - 1)*width(bb)
    dh = 0.5*(s - 1)*height(bb)
    deform(bb, -dw, dw, -dh, dh)
end
(*)(s::Real, bb::BoundingBox) = bb*s

"""
    rotate(bb::BoundingBox, angle, o) -> bbnew

Rotate `bb` around `o` by `angle`, returning the BoundingBox that
encloses the vertices of the rotated box.
"""
function rotate(bb::BoundingBox, angle::Real, p::Point)
    a = rotate(Point(bb.xmin,bb.ymin), angle, p)
    b = rotate(Point(bb.xmax,bb.ymin), angle, p)
    c = rotate(Point(bb.xmin,bb.ymax), angle, p)
    d = rotate(Point(bb.xmax,bb.ymax), angle, p)
    BoundingBox(a, b, c, d)
end

function with_aspect_ratio(bb::BoundingBox, ratio::Real)
    if ratio < aspect_ratio(bb)
        dh = height(bb) - ratio * width(bb)
        return BoundingBox(bb.xmin, bb.xmax, bb.ymin + dh/2, bb.ymax - dh/2)
    else
        dw = width(bb) - height(bb) / ratio
        return BoundingBox(bb.xmin + dw/2, bb.xmax - dw/2, bb.ymin, bb.ymax)
    end
end

"""
    isinside(bb::BoundingBox, p::Point) -> tf::Bool
    isinside(bb::BoundingBox, x, y) -> tf::Bool

Determine whether the point lies within `bb`.
"""
isinside(bb::BoundingBox, x, y) = (bb.xmin <= x <= bb.xmax) && (bb.ymin <= y <= bb.ymax)
isinside(bb::BoundingBox, p::Point) = isinside(bb, p.x, p.y)

# Part 2. Drawing

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

# a graphics output device; can create GraphicsContexts
abstract type GraphicsDevice end

@mustimplement width(gd::GraphicsDevice)
@mustimplement height(gd::GraphicsDevice)
@mustimplement creategc(gd::GraphicsDevice)
xmin(g::GraphicsDevice) = 0
xmax(g::GraphicsDevice) = width(g)
ymin(g::GraphicsDevice) = 0
ymax(g::GraphicsDevice) = height(g)

# an object that can actually be drawn to
abstract type GraphicsContext end

@mustimplement width(gc::GraphicsContext)
@mustimplement height(gc::GraphicsContext)
# getgc() - get a GraphicsContext from something that might be drawable
getgc(gc::GraphicsContext) = gc

# transformations

"""
    inner_canvas(c::GraphicsContext, device::BoundingBox, user::BoundingBox)
    inner_canvas(c::GraphicsContext, x, y, w, h, l, r, t, b)

Create a rectangular drawing area inside `device` (represented in
device-coordinates), giving it user-coordinates `user`. Any drawing
that occurs outside this box is clipped.

`x`, `y`, `w`, and `h` are an alternative parametrization of `device`,
and `l`, `r`, `t`, `b` parametrize `user`.

See also: [`set_coordinates`](@ref).
"""
inner_canvas(c::GraphicsContext, device::BoundingBox, user::BoundingBox) =
    inner_canvas(c,
                 device.xmin, device.ymin, width(device), height(device),
                 user.xmin, user.xmax, user.ymin, user.ymax)

function inner_canvas(c::GraphicsContext, x, y, w, h, l, r, t, b)
    reset_transform(c)
    reset_clip(c)
    rectangle(c, x, y, w, h)
    clip(c)
    _set_coordinates(c, x, y, w, h, l, r, t, b)
end

function set_coordinates(c::GraphicsContext, x, y, w, h, l, r, t, b)
    reset_transform(c)
    _set_coordinates(c, x, y, w, h, l, r, t, b)
end

function _set_coordinates(c::GraphicsContext, x, y, w, h, l, r, t, b)
    if (r-l) != w || (b-t) != h || l != x || t != y
        # note: Cairo assigns integer pixel-space coordinates to the grid
        # points between sample locations, not to the centers of pixels.
        xs = w/(r-l)
        ys = h/(b-t)
        scale(c, xs, ys)

        xcent = (l+r)/2
        ycent = (t+b)/2
        translate(c, -xcent + (w/2 + x)/xs, -ycent + (h/2 + y)/ys)
    end
    c
end

"""
    set_coordinates(c::GraphicsContext, device::BoundingBox, user::BoundingBox)
    set_coordinates(c::GraphicsContext, user::BoundingBox)

Set the device->user coordinate transformation of `c` so that
`device`, expressed in "device coordinates" (pixels), is equivalent to
`user` as expressed in "user coordinates". If `device` is omitted, it
defaults to the full span of `c`,
`BoundingBox(0, width(c), 0, height(c))`.

See also `get_matrix`, `set_matrix`.
"""
set_coordinates(c::GraphicsContext, device::BoundingBox, user::BoundingBox) =
    set_coordinates(c,
               device.xmin, device.ymin, width(device), height(device),
               user.xmin, user.xmax, user.ymin, user.ymax)
set_coordinates(c::GraphicsContext, user::BoundingBox) =
    set_coordinates(c, BoundingBox(0, width(c), 0, height(c)), user)

@mustimplement save(gc::GraphicsContext)
@mustimplement restore(gc::GraphicsContext)
@mustimplement reset_transform(gc::GraphicsContext)
@mustimplement rotate(gc::GraphicsContext, ::Real)
@mustimplement scale(gc::GraphicsContext, ::Real, ::Real)
@mustimplement translate(gc::GraphicsContext, ::Real, ::Real)

user_to_device!(gc::GraphicsContext, c::Vector{Float64}) = c
device_to_user!(gc::GraphicsContext, c::Vector{Float64}) = c
user_to_device_distance!(gc::GraphicsContext, c::Vector{Float64}) = c
device_to_user_distance!(gc::GraphicsContext, c::Vector{Float64}) = c

const d2ubuf = zeros(2)
function device_to_user(gc::GraphicsContext, x::Real, y::Real)
    d2ubuf[1] = x
    d2ubuf[2] = y
    device_to_user!(gc, d2ubuf)
    d2ubuf[1], d2ubuf[2]
end
function user_to_device(gc::GraphicsContext, x::Real, y::Real)
    d2ubuf[1] = x
    d2ubuf[2] = y
    user_to_device!(gc, d2ubuf)
    d2ubuf[1], d2ubuf[2]
end

# drawing and properties

@mustimplement set_line_width(gc::GraphicsContext, ::Real)
@mustimplement set_dash(gc::GraphicsContext, ::Vector{Float64}, ::Real)
@mustimplement set_source_rgb(gc::GraphicsContext, ::Real, ::Real, ::Real)
@mustimplement set_source_rgba(gc::GraphicsContext, ::Real, ::Real, ::Real, ::Real)
@mustimplement set_source(gc::GraphicsContext, src)

# set source color as a Color
function set_source(gc::GraphicsContext, c::Color)
    rgb = convert(RGB, c)
    set_source_rgb(gc, rgb.r, rgb.g, rgb.b)
end

@mustimplement clip(gc::GraphicsContext)
@mustimplement clip_preserve(gc::GraphicsContext)
@mustimplement reset_clip(gc::GraphicsContext)

@mustimplement move_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement line_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement rel_move_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement rel_line_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement arc(gc::GraphicsContext, ::Real, ::Real, ::Real, ::Real, ::Real)
@mustimplement close_path(gc::GraphicsContext)
@mustimplement new_path(gc::GraphicsContext)
@mustimplement new_sub_path(gc::GraphicsContext)

@mustimplement fill(gc::GraphicsContext)
@mustimplement fill_preserve(gc::GraphicsContext)
@mustimplement stroke(gc::GraphicsContext)
@mustimplement stroke_preserve(gc::GraphicsContext)
@mustimplement paint(gc::GraphicsContext)
stroke_transformed(gc::GraphicsContext) = stroke(gc)
stroke_transformed_preserve(gc::GraphicsContext) = stroke_preserve(gc)

# generic path functions

function rectangle(gc::GraphicsContext, x::Real, y::Real, width::Real, height::Real)
    move_to(gc, x, y)
    rel_line_to(gc, width, 0)
    rel_line_to(gc, 0, height)
    rel_line_to(gc, -width, 0)
    close_path(gc)
end
rectangle(gc::GraphicsContext, user::BoundingBox) = rectangle(gc, user.xmin, user.ymin, width(user), height(user))

circle(ctx::GraphicsContext, x::Real, y::Real, r::Real) =
    arc(ctx, x, y, r, 0., 2pi)

function polygon(gc::GraphicsContext, verts::Matrix, idx::Vector)
    move_to(gc, verts[1,idx[1]], verts[2,idx[1]])
    for i=2:length(idx)
        n = idx[i]
        line_to(gc, verts[1,n], verts[2,n])
    end
    close_path(gc)
end

function polygon(self::GraphicsContext, points::AbstractVector)
    move_to(self, points[1].x, points[1].y)
    for i in 2:length(points)
        line_to(self, points[i].x, points[i].y)
    end
    close_path(self)
end

include("deprecations.jl")

end # module