text:
const FLAG_SIMULATION_MIXED_VARY_RANK_SMALL                 = false
const FLAG_SIMULATION_MIXED_VARY_RANK_MEDIUM                = true
const FLAG_SIMULATION_MIXED_VARY_RANK_LARGE                 = false

const FLAG_SIMULATION_MIXED_VARY_MISSING_PERCENTAGE_SMALL   = false
const FLAG_SIMULATION_MIXED_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const FLAG_SIMULATION_MIXED_VARY_MISSING_PERCENTAGE_LARGE   = false

FLAG_SIMULATION_MIXED_VARY_RANK_SMALL ?
  include("simulation_mixed_vary_rank_small.jl")  : nothing

FLAG_SIMULATION_MIXED_VARY_RANK_MEDIUM ?
  include("simulation_mixed_vary_rank_medium.jl")  : nothing


text:
# Default parameter values
const DEFAULT_COARSE_K = 2
const DEFAULT_QUANTIZATION_K = 256
const DEFAULT_QUANTIZATION_M = 1
const DEFAULT_QUANTIZATION_METHOD = :pq
const DEFAULT_COARSE_DISTANCE = Distances.SqEuclidean()
const DEFAULT_COARSE_QUANTIZER = :naive  # :naive or :hnsw
const DEFAULT_QUANTIZATION_DISTANCE = Distances.SqEuclidean()
const DEFAULT_COARSE_MAXITER = 25
const DEFAULT_QUANTIZATION_MAXITER = 25

ext: jl


text:
@testset "collectnamedtuples" begin
    v = [(a = 1, b = 2), (a = 1, b = 3)]
    @test collect_columns(v) == Columns((a = Int[1, 1], b = Int[2, 3]))

    # test inferrability with constant eltype
    itr = [(a = 1, b = 2), (a = 1, b = 2), (a = 1, b = 12)]
    el, st = iterate(itr)
    dest = similar(IndexedTables.arrayof(typeof(el)), 3)
    dest[1] = el
    @inferred IndexedTables.collect_to_columns!(dest, itr, 2, st)

    v = [(a = 1, b = 2), (a = 1.2, b = 3)]
    @test collect_columns(v) == Columns((a = Real[1, 1.2], b = Int[2, 3]))
    @test typeof(collect_columns(v)) == typeof(Columns((a = Real[1, 1.2], b = Int[2, 3])))

    v = [(a = 1, b = 2), (a = 1.2, b = "3")]
    @test collect_columns(v) == Columns((a = Real[1, 1.2], b = Any[2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns((a = Real[1, 1.2], b = Any[2, "3"])))

    v = [(a = 1, b = 2), (a = 1.2, b = 2), (a = 1, b = "3")]
    @test collect_columns(v) == Columns((a = Real[1, 1.2, 1], b = Any[2, 2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns((a = Real[1, 1.2, 1], b = Any[2, 2, "3"])))

    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns((a = [2, 4, 6, 8], b = [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? (a = 1.2, b =i-1) : (a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr_real) == Columns((a = Real[1.2, 4, 6, 8], b = [0, 2, 4, 6]))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns((a = Int[], b = Int[]))

    itr = (i for i in 0:-1)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns((a = Int[], b = Int[]))
end

@testset "collecttuples" begin
    v = [(1, 2), (1, 3)]
    @test collect_columns(v) == Columns((Int[1, 1], Int[2, 3]))
    @inferred collect_columns(v)

    v = [(1, 2), (1.2, 3)]
    @test collect_columns(v) == Columns((Real[1, 1.2], Int[2, 3]))

    v = [(1, 2), (1.2, "3")]
    @test collect_columns(v) == Columns((Real[1, 1.2], Any[2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns((Real[1, 1.2], Any[2, "3"])))

    v = [(1, 2), (1.2, 2), (1, "3")]
    @test collect_columns(v) == Columns((Real[1, 1.2, 1], Any[2, 2, "3"]))
    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(([2, 4, 6, 8], [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? (1.2, i-1) : (i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr_real) == Columns((Real[1.2, 4, 6, 8], [0, 2, 4, 6]))
    @test typeof(collect_columns(tuple_itr_real)) == typeof(Columns((Real[1.2, 4, 6, 8], [0, 2, 4, 6])))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns((Int[], Int[]))

    itr = (i for i in 0:-1)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns((Int[], Int[]))
end

@testset "collectscalars" begin
    v = (i for i in 1:3)
    @test collect_columns(v) == [1,2,3]
    @inferred collect_columns(v)

    v = (i == 1 ? 1.2 : i for i in 1:3)
    @test collect_columns(v) == collect(v)

    itr = Iterators.filter(isodd, 1:100)
    @test collect_columns(itr) == collect(itr)
    real_itr = (i == 1 ? 1.5 : i for i in itr)
    @test collect_columns(real_itr) == collect(real_itr)
    @test eltype(collect_columns(real_itr)) == Real

    #empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = (exp(i) for i in itr)
    @test collect_columns(tuple_itr) == Float64[]

    itr = (i for i in 0:-1)
    tuple_itr = (exp(i) for i in itr)
    @test collect_columns(tuple_itr) == Float64[]

    t = collect_columns((a = i,) for i in (1, missing, 3))
    @test columns(t, 1) isa Vector{Union{Missing, Int}}
    @test isequal(columns(t, 1), [1, missing, 3])
end

@testset "collectpairs" begin
    v = (i=>i+1 for i in 1:3)
    @test collect_columns(v) == columnspair([1,2,3], [2,3,4])
    @test eltype(collect_columns(v)) == Pair{Int, Int}

    v = (i == 1 ? (1.2 => i+1) : (i => i+1) for i in 1:3)
    @test collect_columns(v) == columnspair(Real[1.2,2,3], [2,3,4])
    @test eltype(collect_columns(v)) == Pair{Real, Int}

    v = ((a=i,) => (b="a$i",) for i in 1:3)
    @test collect_columns(v) == columnspair(Columns((a = [1,2,3],)), Columns((b = ["a1","a2","a3"],)))
    @test eltype(collect_columns(v)) == Pair{NamedTuple{(:a,), Tuple{Int}}, NamedTuple{(:b,), Tuple{String}}}

    v = (i == 1 ? (a="1",) => (b="a$i",) : (a=i,) => (b="a$i",) for i in 1:3)
    @test collect_columns(v) == columnspair(Columns((a = ["1",2,3],)), Columns((b = ["a1","a2","a3"],)))
    @test eltype(collect_columns(v)) == Pair{NamedTuple{(:a,), Tuple{Any}}, NamedTuple{(:b,), Tuple{String}}}

    # empty
    v = ((a=i,) => (b="a$i",) for i in 0:-1)
    @test collect_columns(v) == columnspair(Columns((a = Int[],)), Columns((b = String[],)))
    @test eltype(collect_columns(v)) == Pair{NamedTuple{(:a,), Tuple{Int}}, NamedTuple{(:b,), Tuple{String}}}

    v = Iterators.filter(t -> t.first.a == 4, ((a=i,) => (b="a$i",) for i in 1:3))
    @test collect_columns(v) == columnspair(Columns((a = Int[],)), Columns((b = String[],)))
    @test eltype(collect_columns(v)) == Pair{NamedTuple{(:a,), Tuple{Int}}, NamedTuple{(:b,), Tuple{String}}}

    t = table(collect_columns((b = 1,) => (a = i,) for i in (2, missing, 3)))
    @test isequal(t, table((b = [1,1,1], a = [2, missing, 3]), pkey = :b))
end

@testset "collectflattened" begin
    t = [(:a => [1, 2]), (:b => [1, 3])]
    @test collect_columns_flattened(t) == columnspair([:a, :a, :b, :b], [1, 2, 1, 3])
    t = ([(a = 1,), (a = 2,)], [(a = 1.1,), (a = 2.2,)])
    @test collect_columns_flattened(t) == Columns(a = Real[1, 2, 1.1, 2.2])
    @test eltype(collect_columns_flattened(t)) == NamedTuple{(:a,), Tuple{Real}}
    t = [(:a => table(1:2, ["a", "b"])), (:b => table(3:4, ["c", "d"]))]
    @test table(collect_columns_flattened(t)) == table([:a, :a, :b, :b], 1:4, ["a", "b", "c", "d"], pkey = 1)
end


text:
"""
    cross = Cross(c, arm::Int)
A `Drawable` cross passing through the point `c` with arms that are `arm` pixels long.
"""
Cross(c, arm::Int) = Cross(c, -arm:arm)

function draw!(img::AbstractArray{T, 2}, cross::Cross, color::T) where T<:Colorant
    for Δx in cross.range
        drawifinbounds!(img, cross.c.y, cross.c.x + Δx, color)
    end
    for Δy in cross.range
        drawifinbounds!(img, cross.c.y + Δy, cross.c.x, color)
    end
    img
end


text:
module Metrics

"""
    creategrid!(x1, elemtocoord, ξ1)

Create a 1-D grid using `elemtocoord` (see [`brickmesh`](@ref)) using the 1-D
`(-1, 1)` reference coordinates `ξ1`. The element grids are filled using linear
interpolation of the element coordinates.

If `Nq = length(ξ1)` and `nelem = size(elemtocoord, 3)` then the preallocated
array `x1` should be `Nq * nelem == length(x1)`.
"""
function creategrid!(x1, e2c, ξ1)
    (d, nvert, nelem) = size(e2c)
    @assert d == 1
    Nq = length(ξ1)
    x1 = reshape(x1, (Nq, nelem))

    # linear blend
    @inbounds for e in 1:nelem
        for i in 1:Nq
            x1[i, e] =
                ((1 - ξ1[i]) * e2c[1, 1, e] + (1 + ξ1[i]) * e2c[1, 2, e]) / 2
        end
    end
    nothing
end

"""
    creategrid!(x1, x2, elemtocoord, ξ1, ξ2)

Create a 2-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `ξ1` and `ξ2`. The element grids
are filled using bilinear interpolation of the element coordinates.

If `Nq = (length(ξ1), length(ξ2))` and `nelem = size(elemtocoord, 3)` then the
preallocated arrays `x1` and `x2` should be
`prod(Nq) * nelem == size(x1) == size(x2)`.
"""
function creategrid!(x1, x2, e2c, ξ1, ξ2)
    (d, nvert, nelem) = size(e2c)
    @assert d == 2
    Nq = (length(ξ1), length(ξ2))
    x1 = reshape(x1, (Nq..., nelem))
    x2 = reshape(x2, (Nq..., nelem))

    # # bilinear blend of corners
    @inbounds for (f, n) in zip((x1, x2), 1:d)
        for e in 1:nelem, j in 1:Nq[2], i in 1:Nq[1]
            f[i, j, e] =
                (
                    (1 - ξ1[i]) * (1 - ξ2[j]) * e2c[n, 1, e] +
                    (1 + ξ1[i]) * (1 - ξ2[j]) * e2c[n, 2, e] +
                    (1 - ξ1[i]) * (1 + ξ2[j]) * e2c[n, 3, e] +
                    (1 + ξ1[i]) * (1 + ξ2[j]) * e2c[n, 4, e]
                ) / 4
        end
    end
    nothing
end

"""
    creategrid!(x1, x2, x3, elemtocoord, ξ1, ξ2, ξ3)

Create a 3-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `ξ1`. The element grids are filled
using trilinear interpolation of the element coordinates.

If `Nq = (length(ξ1), length(ξ2), length(ξ3))` and
`nelem = size(elemtocoord, 3)` then the preallocated arrays `x1`, `x2`, and `x3`
should be `prod(Nq) * nelem == size(x1) == size(x2) == size(x3)`.
"""
function creategrid!(x1, x2, x3, e2c, ξ1, ξ2, ξ3)
    (d, nvert, nelem) = size(e2c)
    @assert d == 3
    Nq = (length(ξ1), length(ξ2), length(ξ3))
    x1 = reshape(x1, (Nq..., nelem))
    x2 = reshape(x2, (Nq..., nelem))
    x3 = reshape(x3, (Nq..., nelem))

    # trilinear blend of corners
    @inbounds for (f, n) in zip((x1, x2, x3), 1:d)
        for e in 1:nelem, k in 1:Nq[3], j in 1:Nq[2], i in 1:Nq[1]
            f[i, j, k, e] =
                (
                    (1 - ξ1[i]) * (1 - ξ2[j]) * (1 - ξ3[k]) * e2c[n, 1, e] +
                    (1 + ξ1[i]) * (1 - ξ2[j]) * (1 - ξ3[k]) * e2c[n, 2, e] +
                    (1 - ξ1[i]) * (1 + ξ2[j]) * (1 - ξ3[k]) * e2c[n, 3, e] +
                    (1 + ξ1[i]) * (1 + ξ2[j]) * (1 - ξ3[k]) * e2c[n, 4, e] +
                    (1 - ξ1[i]) * (1 - ξ2[j]) * (1 + ξ3[k]) * e2c[n, 5, e] +
                    (1 + ξ1[i]) * (1 - ξ2[j]) * (1 + ξ3[k]) * e2c[n, 6, e] +
                    (1 - ξ1[i]) * (1 + ξ2[j]) * (1 + ξ3[k]) * e2c[n, 7, e] +
                    (1 + ξ1[i]) * (1 + ξ2[j]) * (1 + ξ3[k]) * e2c[n, 8, e]
                ) / 8
        end
    end
    nothing
end

"""
    computemetric!(x1, J, JcV, ξ1x1, sJ, n1, D)

Compute the 1-D metric terms from the element grid arrays `x1`. All the arrays
are preallocated by the user and the (square) derivative matrix `D` should be
consistent with the reference grid `ξ1` used in [`creategrid!`](@ref).

If `Nq = size(D, 1)` and `nelem = div(length(x1), Nq)` then the volume arrays
`x1`, `J`, and `ξ1x1` should all have length `Nq * nelem`.  Similarly, the face
arrays `sJ` and `n1` should be of length `nface * nelem` with `nface = 2`.
"""
function computemetric!(x1, J, JcV, ξ1x1, sJ, n1, D)
    Nq = size(D, 1)
    nelem = div(length(J), Nq)
    x1 = reshape(x1, (Nq, nelem))
    J = reshape(J, (Nq, nelem))
    JcV = reshape(JcV, (Nq, nelem))
    ξ1x1 = reshape(ξ1x1, (Nq, nelem))
    nface = 2
    n1 = reshape(n1, (1, nface, nelem))
    sJ = reshape(sJ, (1, nface, nelem))

    @inbounds for e in 1:nelem
        JcV[:, e] = J[:, e] = D * x1[:, e]
    end
    ξ1x1 .= 1 ./ J

    n1[1, 1, :] .= -sign.(J[1, :])
    n1[1, 2, :] .= sign.(J[Nq, :])
    sJ .= 1
    nothing
end

"""
    computemetric!(x1, x2, J, JcV, ξ1x1, ξ2x1, ξ1x2, ξ2x2, sJ, n1, n2, D1, D2)

Compute the 2-D metric terms from the element grid arrays `x1` and `x2`. All the
arrays are preallocated by the user and the (square) derivative matrice `D1` and
D2 should be consistent with the reference grid `ξ1` and `ξ2` used in
[`creategrid!`](@ref).

If `Nq = (size(D1, 1), size(D2, 1))` and `nelem = div(length(x1), prod(Nq))`
then the volume arrays `x1`, `x2`, `J`, `ξ1x1`, `ξ2x1`, `ξ1x2`, and `ξ2x2`
should all be of size `(Nq..., nelem)`.  Similarly, the face arrays `sJ`, `n1`,
and `n2` should be of size `(maximum(Nq), nface, nelem)` with `nface = 4`
"""
function computemetric!(
    x1,
    x2,
    J,
    JcV,
    ξ1x1,
    ξ2x1,
    ξ1x2,
    ξ2x2,
    sJ,
    n1,
    n2,
    D1,
    D2,
)
    T = eltype(x1)
    Nq = (size(D1, 1), size(D2, 1))
    nelem = div(length(J), prod(Nq))
    x1 = reshape(x1, (Nq..., nelem))
    x2 = reshape(x2, (Nq..., nelem))
    J = reshape(J, (Nq..., nelem))
    JcV = reshape(JcV, (Nq..., nelem))
    ξ1x1 = reshape(ξ1x1, (Nq..., nelem))
    ξ2x1 = reshape(ξ2x1, (Nq..., nelem))
    ξ1x2 = reshape(ξ1x2, (Nq..., nelem))
    ξ2x2 = reshape(ξ2x2, (Nq..., nelem))
    nface = 4
    Nfp = div.(prod(Nq), Nq)
    n1 = reshape(n1, (maximum(Nfp), nface, nelem))
    n2 = reshape(n2, (maximum(Nfp), nface, nelem))
    sJ = reshape(sJ, (maximum(Nfp), nface, nelem))

    for e in 1:nelem
        for j in 1:Nq[2], i in 1:Nq[1]
            x1ξ1 = x1ξ2 = zero(T)
            x2ξ1 = x2ξ2 = zero(T)
            for n in 1:Nq[1]
                x1ξ1 += D1[i, n] * x1[n, j, e]
                x2ξ1 += D1[i, n] * x2[n, j, e]
            end
            for n in 1:Nq[2]
                x1ξ2 += D2[j, n] * x1[i, n, e]
                x2ξ2 += D2[j, n] * x2[i, n, e]
            end
            JcV[i, j, e] = hypot(x1ξ2, x2ξ2)
            J[i, j, e] = x1ξ1 * x2ξ2 - x2ξ1 * x1ξ2
            ξ1x1[i, j, e] = x2ξ2 / J[i, j, e]
            ξ2x1[i, j, e] = -x2ξ1 / J[i, j, e]
            ξ1x2[i, j, e] = -x1ξ2 / J[i, j, e]
            ξ2x2[i, j, e] = x1ξ1 / J[i, j, e]
        end

        for i in 1:maximum(Nfp)
            if i <= Nfp[1]
                n1[i, 1, e] = -J[1, i, e] * ξ1x1[1, i, e]
                n2[i, 1, e] = -J[1, i, e] * ξ1x2[1, i, e]
                n1[i, 2, e] = J[Nq[1], i, e] * ξ1x1[Nq[1], i, e]
                n2[i, 2, e] = J[Nq[1], i, e] * ξ1x2[Nq[1], i, e]
            else
                n1[i, 1:2, e] .= NaN
                n2[i, 1:2, e] .= NaN
            end
            if i <= Nfp[2]
                n1[i, 3, e] = -J[i, 1, e] * ξ2x1[i, 1, e]
                n2[i, 3, e] = -J[i, 1, e] * ξ2x2[i, 1, e]
                n1[i, 4, e] = J[i, Nq[2], e] * ξ2x1[i, Nq[2], e]
                n2[i, 4, e] = J[i, Nq[2], e] * ξ2x2[i, Nq[2], e]
            else
                n1[i, 3:4, e] .= NaN
                n2[i, 3:4, e] .= NaN
            end

            for n in 1:nface
                sJ[i, n, e] = hypot(n1[i, n, e], n2[i, n, e])
                n1[i, n, e] /= sJ[i, n, e]
                n2[i, n, e] /= sJ[i, n, e]
            end
        end
    end

    nothing
end

"""
    computemetric!(x1, x2, x3, J, JcV, ξ1x1, ξ2x1, ξ3x1, ξ1x2, ξ2x2, ξ3x2, ξ1x3,
                   ξ2x3, ξ3x3, sJ, n1, n2, n3, D)

Compute the 3-D metric terms from the element grid arrays `x1`, `x2`, and `x3`.
All the arrays are preallocated by the user and the (square) derivative matrix
`D` should be consistent with the reference grid `ξ1` used in
[`creategrid!`](@ref).

If `Nq = size(D, 1)` and `nelem = div(length(x1), Nq^3)` then the volume arrays
`x1`, `x2`, `x3`, `J`, `ξ1x1`, `ξ2x1`, `ξ3x1`, `ξ1x2`, `ξ2x2`, `ξ3x2`, `ξ1x3`,
`ξ2x3`, and `ξ3x3` should all be of length `Nq^3 * nelem`.  Similarly, the face
arrays `sJ`, `n1`, `n2`, and `n3` should be of size `Nq^2 * nface * nelem` with
`nface = 6`.

The curl invariant formulation of Kopriva (2006), equation 37, is used.

Reference:
 - [Kopriva2006](@cite)
"""
function computemetric!(
    x1,
    x2,
    x3,
    J,
    JcV,
    ξ1x1,
    ξ2x1,
    ξ3x1,
    ξ1x2,
    ξ2x2,
    ξ3x2,
    ξ1x3,
    ξ2x3,
    ξ3x3,
    sJ,
    n1,
    n2,
    n3,
    D1,
    D2,
    D3,
)
    T = eltype(x1)

    Nq = (size(D1, 1), size(D2, 1), size(D3, 1))
    Np = prod(Nq)
    Nfp = div.(Np, Nq)
    nelem = div(length(J), Np)

    x1 = reshape(x1, (Nq..., nelem))
    x2 = reshape(x2, (Nq..., nelem))
    x3 = reshape(x3, (Nq..., nelem))
    J = reshape(J, (Nq..., nelem))
    JcV = reshape(JcV, (Nq..., nelem))
    ξ1x1 = reshape(ξ1x1, (Nq..., nelem))
    ξ2x1 = reshape(ξ2x1, (Nq..., nelem))
    ξ3x1 = reshape(ξ3x1, (Nq..., nelem))
    ξ1x2 = reshape(ξ1x2, (Nq..., nelem))
    ξ2x2 = reshape(ξ2x2, (Nq..., nelem))
    ξ3x2 = reshape(ξ3x2, (Nq..., nelem))
    ξ1x3 = reshape(ξ1x3, (Nq..., nelem))
    ξ2x3 = reshape(ξ2x3, (Nq..., nelem))
    ξ3x3 = reshape(ξ3x3, (Nq..., nelem))

    nface = 6
    n1 = reshape(n1, maximum(Nfp), nface, nelem)
    n2 = reshape(n2, maximum(Nfp), nface, nelem)
    n3 = reshape(n3, maximum(Nfp), nface, nelem)
    sJ = reshape(sJ, maximum(Nfp), nface, nelem)

    JI2 = similar(J, Nq...)
    (yzr, yzs, yzt) = (similar(JI2), similar(JI2), similar(JI2))
    (zxr, zxs, zxt) = (similar(JI2), similar(JI2), similar(JI2))
    (xyr, xys, xyt) = (similar(JI2), similar(JI2), similar(JI2))

    ξ1x1 .= zero(T)
    ξ2x1 .= zero(T)
    ξ3x1 .= zero(T)
    ξ1x2 .= zero(T)
    ξ2x2 .= zero(T)
    ξ3x2 .= zero(T)
    ξ1x3 .= zero(T)
    ξ2x3 .= zero(T)
    ξ3x3 .= zero(T)

    fill!(n1, NaN)
    fill!(n2, NaN)
    fill!(n3, NaN)
    fill!(sJ, NaN)

    @inbounds for e in 1:nelem
        for k in 1:Nq[3], j in 1:Nq[2], i in 1:Nq[1]
            x1ξ1 = x1ξ2 = x1ξ3 = zero(T)
            x2ξ1 = x2ξ2 = x2ξ3 = zero(T)
            x3ξ1 = x3ξ2 = x3ξ3 = zero(T)
            for n in 1:Nq[1]
                x1ξ1 += D1[i, n] * x1[n, j, k, e]
                x2ξ1 += D1[i, n] * x2[n, j, k, e]
                x3ξ1 += D1[i, n] * x3[n, j, k, e]
            end
            for n in 1:Nq[2]
                x1ξ2 += D2[j, n] * x1[i, n, k, e]
                x2ξ2 += D2[j, n] * x2[i, n, k, e]
                x3ξ2 += D2[j, n] * x3[i, n, k, e]
            end
            for n in 1:Nq[3]
                x1ξ3 += D3[k, n] * x1[i, j, n, e]
                x2ξ3 += D3[k, n] * x2[i, j, n, e]
                x3ξ3 += D3[k, n] * x3[i, j, n, e]
            end
            JcV[i, j, k, e] = hypot(x1ξ3, x2ξ3, x3ξ3)
            J[i, j, k, e] = (
                x1ξ1 * (x2ξ2 * x3ξ3 - x3ξ2 * x2ξ3) +
                x2ξ1 * (x3ξ2 * x1ξ3 - x1ξ2 * x3ξ3) +
                x3ξ1 * (x1ξ2 * x2ξ3 - x2ξ2 * x1ξ3)
            )

            JI2[i, j, k] = 1 / (2 * J[i, j, k, e])

            yzr[i, j, k] = x2[i, j, k, e] * x3ξ1 - x3[i, j, k, e] * x2ξ1
            yzs[i, j, k] = x2[i, j, k, e] * x3ξ2 - x3[i, j, k, e] * x2ξ2
            yzt[i, j, k] = x2[i, j, k, e] * x3ξ3 - x3[i, j, k, e] * x2ξ3
            zxr[i, j, k] = x3[i, j, k, e] * x1ξ1 - x1[i, j, k, e] * x3ξ1
            zxs[i, j, k] = x3[i, j, k, e] * x1ξ2 - x1[i, j, k, e] * x3ξ2
            zxt[i, j, k] = x3[i, j, k, e] * x1ξ3 - x1[i, j, k, e] * x3ξ3
            xyr[i, j, k] = x1[i, j, k, e] * x2ξ1 - x2[i, j, k, e] * x1ξ1
            xys[i, j, k] = x1[i, j, k, e] * x2ξ2 - x2[i, j, k, e] * x1ξ2
            xyt[i, j, k] = x1[i, j, k, e] * x2ξ3 - x2[i, j, k, e] * x1ξ3
        end

        for k in 1:Nq[3], j in 1:Nq[2], i in 1:Nq[1]
            for n in 1:Nq[1]
                ξ2x1[i, j, k, e] -= D1[i, n] * yzt[n, j, k]
                ξ3x1[i, j, k, e] += D1[i, n] * yzs[n, j, k]
                ξ2x2[i, j, k, e] -= D1[i, n] * zxt[n, j, k]
                ξ3x2[i, j, k, e] += D1[i, n] * zxs[n, j, k]
                ξ2x3[i, j, k, e] -= D1[i, n] * xyt[n, j, k]
                ξ3x3[i, j, k, e] += D1[i, n] * xys[n, j, k]
            end
            for n in 1:Nq[2]
                ξ1x1[i, j, k, e] += D2[j, n] * yzt[i, n, k]
                ξ3x1[i, j, k, e] -= D2[j, n] * yzr[i, n, k]
                ξ1x2[i, j, k, e] += D2[j, n] * zxt[i, n, k]
                ξ3x2[i, j, k, e] -= D2[j, n] * zxr[i, n, k]
                ξ1x3[i, j, k, e] += D2[j, n] * xyt[i, n, k]
                ξ3x3[i, j, k, e] -= D2[j, n] * xyr[i, n, k]
            end
            for n in 1:Nq[3]
                ξ1x1[i, j, k, e] -= D3[k, n] * yzs[i, j, n]
                ξ2x1[i, j, k, e] += D3[k, n] * yzr[i, j, n]
                ξ1x2[i, j, k, e] -= D3[k, n] * zxs[i, j, n]
                ξ2x2[i, j, k, e] += D3[k, n] * zxr[i, j, n]
                ξ1x3[i, j, k, e] -= D3[k, n] * xys[i, j, n]
                ξ2x3[i, j, k, e] += D3[k, n] * xyr[i, j, n]
            end
            ξ1x1[i, j, k, e] *= JI2[i, j, k]
            ξ2x1[i, j, k, e] *= JI2[i, j, k]
            ξ3x1[i, j, k, e] *= JI2[i, j, k]
            ξ1x2[i, j, k, e] *= JI2[i, j, k]
            ξ2x2[i, j, k, e] *= JI2[i, j, k]
            ξ3x2[i, j, k, e] *= JI2[i, j, k]
            ξ1x3[i, j, k, e] *= JI2[i, j, k]
            ξ2x3[i, j, k, e] *= JI2[i, j, k]
            ξ3x3[i, j, k, e] *= JI2[i, j, k]
        end

        # faces 1 & 2
        for k in 1:Nq[3], j in 1:Nq[2]
            n = j + (k - 1) * Nq[2]
            n1[n, 1, e] = -J[1, j, k, e] * ξ1x1[1, j, k, e]
            n2[n, 1, e] = -J[1, j, k, e] * ξ1x2[1, j, k, e]
            n3[n, 1, e] = -J[1, j, k, e] * ξ1x3[1, j, k, e]
            n1[n, 2, e] = J[Nq[1], j, k, e] * ξ1x1[Nq[1], j, k, e]
            n2[n, 2, e] = J[Nq[1], j, k, e] * ξ1x2[Nq[1], j, k, e]
            n3[n, 2, e] = J[Nq[1], j, k, e] * ξ1x3[Nq[1], j, k, e]
            for f in 1:2
                sJ[n, f, e] = hypot(n1[n, f, e], n2[n, f, e], n3[n, f, e])
                n1[n, f, e] /= sJ[n, f, e]
                n2[n, f, e] /= sJ[n, f, e]
                n3[n, f, e] /= sJ[n, f, e]
            end
        end
        # faces 3 & 4
        for k in 1:Nq[3], i in 1:Nq[1]
            n = i + (k - 1) * Nq[1]
            n1[n, 3, e] = -J[i, 1, k, e] * ξ2x1[i, 1, k, e]
            n2[n, 3, e] = -J[i, 1, k, e] * ξ2x2[i, 1, k, e]
            n3[n, 3, e] = -J[i, 1, k, e] * ξ2x3[i, 1, k, e]
            n1[n, 4, e] = J[i, Nq[2], k, e] * ξ2x1[i, Nq[2], k, e]
            n2[n, 4, e] = J[i, Nq[2], k, e] * ξ2x2[i, Nq[2], k, e]
            n3[n, 4, e] = J[i, Nq[2], k, e] * ξ2x3[i, Nq[2], k, e]
            for f in 3:4
                sJ[n, f, e] = hypot(n1[n, f, e], n2[n, f, e], n3[n, f, e])
                n1[n, f, e] /= sJ[n, f, e]
                n2[n, f, e] /= sJ[n, f, e]
                n3[n, f, e] /= sJ[n, f, e]
            end
        end
        # faces 5 & 6
        for j in 1:Nq[2], i in 1:Nq[1]
            n = i + (j - 1) * Nq[1]
            n1[n, 5, e] = -J[i, j, 1, e] * ξ3x1[i, j, 1, e]
            n2[n, 5, e] = -J[i, j, 1, e] * ξ3x2[i, j, 1, e]
            n3[n, 5, e] = -J[i, j, 1, e] * ξ3x3[i, j, 1, e]
            n1[n, 6, e] = J[i, j, Nq[3], e] * ξ3x1[i, j, Nq[3], e]
            n2[n, 6, e] = J[i, j, Nq[3], e] * ξ3x2[i, j, Nq[3], e]
            n3[n, 6, e] = J[i, j, Nq[3], e] * ξ3x3[i, j, Nq[3], e]
            for f in 5:6
                sJ[n, f, e] = hypot(n1[n, f, e], n2[n, f, e], n3[n, f, e])
                n1[n, f, e] /= sJ[n, f, e]
                n2[n, f, e] /= sJ[n, f, e]
                n3[n, f, e] /= sJ[n, f, e]
            end
        end
    end

    nothing
end

"""
  computemetric(x1::AbstractArray{T, 2}, D::AbstractMatrix{T}) where T

Compute the 1-D metric terms from the element grid array `x1` using the
derivative matrix `D`. The derivative matrix `D` should be consistent with the
reference grid `ξ1` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

- `J` the Jacobian determinant
- `ξ1x1` derivative ∂ξ1 / ∂x1'
- `sJ` the surface Jacobian
- 'n1` outward pointing unit normal in \$x1\$-direction
"""
function computemetric(x1::AbstractArray{T, 2}, D::AbstractMatrix{T}) where {T}

    Nq = size(D, 1)
    nelem = size(x1, 2)
    nface = 2

    J = similar(x1)
    ξ1x1 = similar(x1)

    sJ = Array{T, 3}(undef, 1, nface, nelem)
    n1 = Array{T, 3}(undef, 1, nface, nelem)

    computemetric!(x1, J, JcV, ξ1x1, sJ, n1, D)

    (J = J, JcV = JcV, ξ1x1 = ξ1x1, sJ = sJ, n1 = n1)
end

"""
  computemetric(x1::AbstractArray{T, 3}, x2::AbstractArray{T, 3},
                D::AbstractMatrix{T}) where T

Compute the 2-D metric terms from the element grid arrays `x1` and `x2` using
the derivative matrix `D`. The derivative matrix `D` should be consistent with
the reference grid `ξ1` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

- `J` the Jacobian determinant
- `ξ1x1` derivative ∂ξ1 / ∂x1'
- `ξ2x1` derivative ∂ξ2 / ∂x1'
 - `ξ1x2` derivative ∂ξ1 / ∂x2'
 - `ξ2x2` derivative ∂ξ2 / ∂x2'
 - `sJ` the surface Jacobian
 - 'n1` outward pointing unit normal in \$x1\$-direction
 - 'n2` outward pointing unit normal in \$x2\$-direction
"""
function computemetric(
    x1::AbstractArray{T, 3},
    x2::AbstractArray{T, 3},
    D::AbstractMatrix{T},
) where {T}
    @assert size(x1) == size(x2)
    Nq = size(D, 1)
    nelem = size(x1, 3)
    nface = 4

    J = similar(x1)
    JcV = similar(x1)
    ξ1x1 = similar(x1)
    ξ2x1 = similar(x1)
    ξ1x2 = similar(x1)
    ξ2x2 = similar(x1)

    sJ = Array{T, 3}(undef, Nq, nface, nelem)
    n1 = Array{T, 3}(undef, Nq, nface, nelem)
    n2 = Array{T, 3}(undef, Nq, nface, nelem)

    computemetric!(x1, x2, J, JcV, ξ1x1, ξ2x1, ξ1x2, ξ2x2, sJ, n1, n2, D)

    (
        J = J,
        JcV = JcV,
        ξ1x1 = ξ1x1,
        ξ2x1 = ξ2x1,
        ξ1x2 = ξ1x2,
        ξ2x2 = ξ2x2,
        sJ = sJ,
        n1 = n1,
        n2 = n2,
    )
end

"""
    computemetric(x1::AbstractArray{T, 3}, x2::AbstractArray{T, 3},
                  D::AbstractMatrix{T}) where T

Compute the 3-D metric terms from the element grid arrays `x1`, `x2`, and `x3`
using the derivative matrix `D`. The derivative matrix `D` should be consistent
with the reference grid `ξ1` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

 - `J` the Jacobian determinant
 - `ξ1x1` derivative ∂ξ1 / ∂x1'
 - `ξ2x1` derivative ∂ξ2 / ∂x1'
 - `ξ3x1` derivative ∂ξ3 / ∂x1'
 - `ξ1x2` derivative ∂ξ1 / ∂x2'
 - `ξ2x2` derivative ∂ξ2 / ∂x2'
 - `ξ3x2` derivative ∂ξ3 / ∂x2'
 - `ξ1x3` derivative ∂ξ1 / ∂x3'
 - `ξ2x3` derivative ∂ξ2 / ∂x3'
 - `ξ3x3` derivative ∂ξ3 / ∂x3'
 - `sJ` the surface Jacobian
 - 'n1` outward pointing unit normal in \$x1\$-direction
 - 'n2` outward pointing unit normal in \$x2\$-direction
 - 'n3` outward pointing unit normal in \$x3\$-direction

!!! note

   The storage of the volume terms and surface terms from this function are
   slightly different. The volume terms used Cartesian indexing whereas the
   surface terms use linear indexing.
"""
function computemetric(
    x1::AbstractArray{T, 4},
    x2::AbstractArray{T, 4},
    x3::AbstractArray{T, 4},
    D::AbstractMatrix{T},
) where {T}

    @assert size(x1) == size(x2) == size(x3)
    Nq = size(D, 1)
    nelem = size(x1, 4)
    nface = 6

    J = similar(x1)
    JcV = similar(x1)
    ξ1x1 = similar(x1)
    ξ2x1 = similar(x1)
    ξ3x1 = similar(x1)
    ξ1x2 = similar(x1)
    ξ2x2 = similar(x1)
    ξ3x2 = similar(x1)
    ξ1x3 = similar(x1)
    ξ2x3 = similar(x1)
    ξ3x3 = similar(x1)

    sJ = Array{T, 3}(undef, Nq^2, nface, nelem)
    n1 = Array{T, 3}(undef, Nq^2, nface, nelem)
    n2 = Array{T, 3}(undef, Nq^2, nface, nelem)
    n3 = Array{T, 3}(undef, Nq^2, nface, nelem)

    computemetric!(
        x1,
        x2,
        x3,
        J,
        JcV,
        ξ1x1,
        ξ2x1,
        ξ3x1,
        ξ1x2,
        ξ2x2,
        ξ3x2,
        ξ1x3,
        ξ2x3,
        ξ3x3,
        sJ,
        n1,
        n2,
        n3,
        D,
    )

    (
        J = J,
        JcV = JcV,
        ξ1x1 = ξ1x1,
        ξ2x1 = ξ2x1,
        ξ3x1 = ξ3x1,
        ξ1x2 = ξ1x2,
        ξ2x2 = ξ2x2,
        ξ3x2 = ξ3x2,
        ξ1x3 = ξ1x3,
        ξ2x3 = ξ2x3,
        ξ3x3 = ξ3x3,
        sJ = sJ,
        n1 = n1,
        n2 = n2,
        n3 = n3,
    )
end

end # module

max_forks_repo_name: jm-c/CLIMA


text:
using Game2048Core
using Test

@testset "Game2048Core.jl" begin
    bitboard = initbboard()

    move(bitboard, up)
    move(bitboard, down)
    move(bitboard, left)
    move(bitboard, down)

    simulate_bb(bitboard)

    simulate_bb()
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Makie.backend_display, (CairoBackend, Scene))
    activate!()
    f, ax1, pl = scatter(1:4)
    f, ax2, pl = lines(1:4)
    Makie.colorbuffer(ax1.scene)
    Makie.colorbuffer(ax2.scene)
    return
end


text:
using ACTRSimulators, Test, ACTRModels, Random, DataFrames
import ACTRSimulators: start!, press_key!
Random.seed!(8985)
include("task.jl")

scheduler = ACTRScheduler(;model_trace=true, store=true)
task = SimpleTask(;scheduler)
procedural = Procedural()
T = vo_to_chunk() |> typeof
visual_location = VisualLocation(buffer=T[])
visual = Visual(buffer=T[])
visicon = VisualObject[]
motor = Motor()
memory = [Chunk(;animal=:dog), Chunk(;animal=:cat)]
declarative = Declarative(;memory)
actr = ACTR(;scheduler, procedural, visual_location, visual, motor, declarative, visicon)

function can_attend(actr)
    c1(actr) = !isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_respond(actr)
    c1(actr) = !isempty(actr.visual.buffer)
    c2(actr) = !actr.motor.state.busy
    c3(actr) = !actr.imaginal.state.busy
    return c1, c2, c3
end

function attend_action(actr, task)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk)
    return nothing
end

function motor_action(actr, task)
    buffer = actr.visual.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual)
    encoding!(actr, chunk)

    key = chunk.slots.text
    responding!(actr, task, key)
    return nothing
end

rule1 = Rule(;conditions=can_attend, action=attend_action, actr, task, name="Attend")
push!(procedural.rules, rule1)

rule2 = Rule(;conditions=can_respond, action=motor_action, actr, task, name="Respond")
push!(procedural.rules, rule2)

run!(actr, task)

@test isempty(task.screen)
@test !isempty(actr.imaginal.buffer)

observed = map(x->x.description, scheduler.complete_events)
expected = [
    "Starting",
    "Present Stimulus",
    "Selected Attend",
    "Attend",
    "Selected Respond",
    "Respond",
    "Create New Chunk"
]
@test expected == observed


text:
cd(@__DIR__)
using Pkg
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.activate(@__DIR__)
CI && Pkg.instantiate()
ENV["GKSwstype"] = "100" # allow local builds without output
using DelayEmbeddings
using TransferEntropy
using Documenter
using DocumenterTools: Themes
using CausalityTools
using DynamicalSystems
using HypothesisTests
using Distributions

# %% JuliaDynamics theme.
# download the themes
using DocumenterTools: Themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

# %% Build docs
cd(@__DIR__)
ENV["JULIA_DEBUG"] = "Documenter"

PAGES = [
    "Overview" => "index.md",
    "surrogate.md",

    "Distance based" => [
        "joint_distance_distribution.md",
        "s_measure.md",
        "cross_mapping.md",
        "pairwise_asymmetric_inference.md"
    ],
    "Information/entropy based" => [
        "mutualinfo.md",
        "conditional_mutualinfo.md",
        "TransferEntropy.md",
        "predictive_asymmetry.md",
        "generalized_entropy.md",
        "info_estimators.md",
    ],

    "example_systems.md",
    "Utilities" => [
        "invariant_measure.md",
        "dataset.md",
    ],
 ]

makedocs(
    modules = [CausalityTools, Entropies, TransferEntropy, DelayEmbeddings],
    format = Documenter.HTML(
        prettyurls = CI,
        assets = [
            asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
        ],
        ),
    sitename = "CausalityTools.jl",
    authors = "Kristian Agasøster Haaga, Tor Einar Møller, George Datseris",
    pages = PAGES
)

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/CausalityTools.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
allequal(x) = true
allequal(x, y, z...) = x == y && allequal(y, z...)
function Base.map!(f, y::GPUArray, xs::GPUArray...)
    @assert allequal(size.((y, xs...))...)
    return y .= f.(xs...)
end
function Base.map(f, y::GPUArray, xs::GPUArray...)
    @assert allequal(size.((y, xs...))...)
    return f.(y, xs...)
end

# Break ambiguities with base
Base.map!(f, y::GPUArray) =
    invoke(map!, Tuple{Any,GPUArray,Vararg{GPUArray}}, f, y)
Base.map!(f, y::GPUArray, x::GPUArray) =
    invoke(map!, Tuple{Any,GPUArray, Vararg{GPUArray}}, f, y, x)
Base.map!(f, y::GPUArray, x1::GPUArray, x2::GPUArray) =
    invoke(map!, Tuple{Any,GPUArray, Vararg{GPUArray}}, f, y, x1, x2)

# TODO find out why this segfaults julia without stack trace on AMD
# produces wrong results on Titan X and passes on GTX 950..........

# @generated function nindex(i::T, ls::NTuple{N}) where {T, N}
#     quote
#         Base.@_inline_meta
#         $(foldr(:($T(0), $T(0)), T(1):T(N)) do n, els
#             :(i ≤ ls[$n] ? ($T($n), i) : (i -= $T(ls[$n]); $els))
#         end)
#     end
# end
# function catindex(dim, I::NTuple{N, T}, shapes) where {T, N}
#     xi = nindex(I[dim], map(s-> s[dim], shapes))
#     x = xi[1]; i = xi[2]
#     x, ntuple(n -> n == dim ? i : I[n], Val{N})
# end
#
# function _cat(dim, dest, xs...)
#     gpu_call(dest, (Int(dim), dest, xs)) do state, dim, dest, xs
#         I = @cartesianidx dest state
#         nI = catindex(dim, I, size.(xs))
#         n = nI[1]; I′ = nI[2]
#         @inbounds dest[I...] = xs[n][I′...]
#         return
#     end
#     return dest
# end
#
# function Base.cat_t(dims::Integer, T::Type, x::GPUArray, xs::GPUArray...)
#     catdims = Base.dims2cat(dims)
#     shape = Base.cat_shape(catdims, (), size.((x, xs...))...)
#     dest = Base.cat_similar(x, T, shape)
#     _cat(dims, dest, x, xs...)
# end
#
# Base.vcat(xs::GPUArray...) = cat(1, xs...)
# Base.hcat(xs::GPUArray...) = cat(2, xs...)

# Base functions that are sadly not fit for the the GPU yet (they only work for Int64)
Base.@pure @inline function gpu_ind2sub(A::AbstractArray, ind::T) where T
    _ind2sub(size(A), ind - T(1))
end
Base.@pure @inline function gpu_ind2sub(dims::NTuple{N}, ind::T) where {N, T}
    _ind2sub(NTuple{N, T}(dims), ind - T(1))
end
Base.@pure @inline _ind2sub(::Tuple{}, ind::T) where {T} = (ind + T(1),)
Base.@pure @inline function _ind2sub(indslast::NTuple{1}, ind::T) where T
    ((ind + T(1)),)
end
Base.@pure @inline function _ind2sub(inds, ind::T) where T
    r1 = inds[1]
    indnext = div(ind, r1)
    f = T(1); l = r1
    (ind-l*indnext+f, _ind2sub(Base.tail(inds), indnext)...)
end

Base.@pure function gpu_sub2ind(dims::NTuple{N}, I::NTuple{N2, T}) where {N, N2, T}
    Base.@_inline_meta
    _sub2ind(NTuple{N, T}(dims), T(1), T(1), I...)
end
_sub2ind(x, L, ind) = ind
function _sub2ind(::Tuple{}, L, ind, i::T, I::T...) where T
    Base.@_inline_meta
    ind + (i - T(1)) * L
end
function _sub2ind(inds, L, ind, i::IT, I::IT...) where IT
    Base.@_inline_meta
    r1 = inds[1]
    _sub2ind(Base.tail(inds), L * r1, ind + (i - IT(1)) * L, I...)
end

# This is pretty ugly, but I feel bad to add those to device arrays, since
# we're never bound checking... So getindex(a::GPUVector, 10, 10) would silently go unnoticed
# we need this here for easier implementation of repeat
@inline Base.@propagate_inbounds getidx_2d1d(x::AbstractVector, i, j) = x[i]
@inline Base.@propagate_inbounds getidx_2d1d(x::AbstractMatrix, i, j) = x[i, j]

function Base.repeat(a::GPUVecOrMat, m::Int, n::Int = 1)
    o, p = size(a, 1), size(a, 2)
    b = similar(a, o*m, p*n)
    gpu_call(a, (b, a, o, p, m, n), n) do state, b, a, o, p, m, n
        j = linear_index(state)
        j > n && return
        d = (j - 1) * p + 1
        @inbounds for i in 1:m
            c = (i - 1) * o + 1
            for r in 1:p
                for k in 1:o
                    b[k - 1 + c, r - 1 + d] = getidx_2d1d(a, k, r)
                end
            end
        end
        return
    end
    return b
end

function Base.repeat(a::GPUVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    gpu_call(a, (b, a, o, m), m) do state, b, a, o, m
        i = linear_index(state)
        i > m && return
        c = (i - 1)*o + 1
        @inbounds for i in 1:o
            b[c + i - 1] = a[i]
        end
        return
    end
    return b
end


text:
@testset "Pairwise" begin
  Γ = pairwise(GaussianVariogram(), Matrix(1.0I, 3, 3))
  @test eltype(Γ) == Float64
  @test issymmetric(Γ)

  Γ_f = pairwise(GaussianVariogram(range=1f0, sill=1f0, nugget=0f0), Matrix(1.0f0I, 3, 3))
  @test eltype(Γ_f) == Float32
  @test issymmetric(Γ_f)

  grid = RegularGrid{Float64}(10, 10)
  Γ = pairwise(GaussianVariogram(), grid, 1:5)
  @test size(Γ) == (5, 5)
  @test issymmetric(Γ)
  Γ = pairwise(GaussianVariogram(), grid, 1:3, 7:10)
  @test size(Γ) == (3, 4)
  @test all(Γ .> 0)
end


text:
"""
    DubinsCar

A simple unicyle model with state vector  `[x,y,θ]` and control vector `[v,ω]`.

# Constructor
    DubinsCar(; [radius=0.175])

where `radius` is the radius of the car, used only for plotting and obstacle avoidance constraints.
"""
DubinsCar
@autodiff struct DubinsCar <: ContinuousDynamics
    radius::Float64  # [m] radius of a Roomba
    function DubinsCar(;radius=0.175)
        new(radius)
    end
end
RobotDynamics.state_dim(::DubinsCar) = 3
RobotDynamics.control_dim(::DubinsCar) = 2

function dynamics(::DubinsCar,x,u)
    ẋ = @SVector [u[1]*cos(x[3]),
                  u[1]*sin(x[3]),
                  u[2]]
end
function dynamics!(::DubinsCar,ẋ,x,u)
    ẋ[1] = u[1]*cos(x[3])
    ẋ[2] = u[1]*sin(x[3])
    ẋ[3] = u[2]
    return nothing
end

function jacobian!(::DubinsCar, J,xdot,x,u)
    sθ,cθ = sincos(x[3])
    J .= 0
    J[1,3] = -u[1]*sθ
    J[1,4] = cθ
    J[2,3] = u[1]*cθ
    J[2,4] = sθ
    J[3,5] = 1.0
    return nothing
end

Base.position(::DubinsCar, x) = @SVector [x[1], x[2], 0.0]
orientation(::DubinsCar, x) = expm(x[3]*@SVector [0,0,1.])

"""
    BicycleModel <: ContinuousDynamics

Kinematic model of a car with front-wheel steering. Assumes no skid and models wheels on the same
axel as a single wheel. The state is defined as `[x,y,θ,δ]` where `x` and `y` are the position,
`θ` is the orientation (yaw angle), and `δ` is the steering angle. The controls are `[v,ϕ]` where
`v` is the forward velocity and `ϕ` is the steering angle rate.

The reference location on the car can be specified using the `ref` field, which is either at
the center of the rear wheel `:rear`, the front wheel `:front`, or at the center of gravity `:cg` (default).

# Constructor
    BicycleModel(;kwargs...)

with keyword arguments:
* `ref`: location of the reference position on the vehicle. Valid values: `[:rear, :front, :cg]`
* `L`: distance between the front and rear wheels
* `lr`: distance from the center of gravity to the center of the rear wheel
"""
@autodiff struct BicycleModel <: ContinuousDynamics
    ref::Symbol
    L::Float64    # distance between wheels  (m)
    lr::Float64   # distance to rear wheels  (m)
    function BicycleModel(;ref::Symbol=:cg, L::Real=2.7, lr::Real=1.5)
        @assert ref ∈ (:cg, :front, :rear)
        @assert L > 0 "Vehicle length ($L) must be greater than 0"
        @assert L > lr "Vehicle length ($L) must be greater than the distance to the rear wheel ($lr)"
        new(ref, L, lr)
    end
end

fun = :foo
args = [:(model::BicycleModel), :x, :u]

body = quote
    v = u[1]  # longitudinal velocity (cmd)
    ϕ = u[2]  # steering angle rate   (cmd)
    θ = x[3]  # yaw
    δ = x[4]  # steering angle
    if model.ref == :cg
        β = atan(model.lr * δ, model.L)
        s,c = sincos(θ + β)
        ω = v*cos(β)*tan(δ) / model.L
    elseif model.ref == :rear
        s,c = sincos(θ)
        ω = v*tan(δ) / model.L
    elseif model.ref == :front
        s,c = sincos(θ + δ)
        ω = v*sin(δ) / model.L
    end
    ẋ = v*c
    ẏ = v*s
end
@eval function dynamics(model::BicycleModel, x, u)
    $body
    return SA[ẋ, ẏ, ω, ϕ]
end
@eval function dynamics!(model::BicycleModel, xdot, x, u)
    $body
    xdot[1] = ẋ
    xdot[2] = ẏ
    xdot[3] = ω
    xdot[4] = ϕ
    return nothing
end

state_dim(::BicycleModel) = 4
control_dim(::BicycleModel) = 2

# using Plots
# function Plots.plot(model::DubinsCar, Z::Traj; kwargs...)
#     x = [z.z[1] for z in Z]
#     y = [z.z[2] for z in Z]
#     Plots.plot(x, y, xlabel="x position", ylabel="y position", markershape=:circle;
#         kwargs...)
# end
# function Plots.plot!(model::DubinsCar, Z::Traj; kwargs...)
#     x = [z.z[1] for z in Z]
#     y = [z.z[2] for z in Z]
#     Plots.plot!(x, y, xlabel="x position", ylabel="y position", markershape=:circle;
#         kwargs...)
# end


text:
module Meshing

using GeometryTypes,
      StaticArrays

"""
    _DEFAULT_SAMPLES = (24,24,24)

Global default sampling count for functions.
"""
const _DEFAULT_SAMPLES = (24,24,24)

include("algorithmtypes.jl")
include("geometrytypes_api.jl")
include("common.jl")
include("marching_tetrahedra.jl")
include("marching_cubes.jl")
include("surface_nets.jl")

export isosurface,
       MarchingCubes,
       MarchingTetrahedra,
       NaiveSurfaceNets

end # module


text:
"""
    mexs = parse_source(filename::AbstractString, mod::Module)

Parse the source `filename`, returning a [`ModuleExprsSigs`](@ref) `mexs`.
`mod` is the "parent" module for the file (i.e., the one that `include`d the file);
if `filename` defines more module(s) then these will all have separate entries in `mexs`.

If parsing `filename` fails, `nothing` is returned.
"""
parse_source(filename::AbstractString, mod::Module) =
    parse_source!(ModuleExprsSigs(mod), filename, mod)

"""
    parse_source!(mexs::ModuleExprsSigs, filename, mod::Module)

Top-level parsing of `filename` as included into module
`mod`. Successfully-parsed expressions will be added to `mexs`. Returns
`mexs` if parsing finished successfully, otherwise `nothing` is returned.

See also [`Revise.parse_source`](@ref).
"""
function parse_source!(mod_exprs_sigs::ModuleExprsSigs, filename::AbstractString, mod::Module)
    if !isfile(filename)
        @warn "$filename is not a file, omitting from revision tracking"
        return nothing
    end
    parse_source!(mod_exprs_sigs, read(filename, String), filename, mod)
end

"""
    success = parse_source!(mod_exprs_sigs::ModuleExprsSigs, src::AbstractString, filename::AbstractString, mod::Module)

Parse a string `src` obtained by reading `file` as a single
string. `pos` is the 1-based byte offset from which to begin parsing `src`.

See also [`Revise.parse_source`](@ref).
"""
function parse_source!(mod_exprs_sigs::ModuleExprsSigs, src::AbstractString, filename::AbstractString, mod::Module)
    ex = Base.parse_input_line(src; filename=filename)
    ex === nothing && return mod_exprs_sigs
    if isexpr(ex, :error) || isexpr(ex, :incomplete)
        throw(LoadError(filename, 0, ex.args[1]))  # would be nice to get the line number
    end
    modexs, docexprs = Tuple{Module,Expr}[], DocExprs()
    JuliaInterpreter.split_expressions!(modexs, docexprs, mod, ex; extract_docexprs=true)
    for (mod, ex) in modexs
        exprs_sigs = get(mod_exprs_sigs, mod, nothing)
        if exprs_sigs === nothing
            mod_exprs_sigs[mod] = exprs_sigs = ExprsSigs()
        end
        exprs_sigs[unwrap(ex)] = nothing
        # exprs_sigs[ex] = nothing
    end
    for (mod, docexs) in docexprs
        exprs_sigs = get(mod_exprs_sigs, mod, nothing)
        for dex in docexs
            if length(dex.args) < 5
                push!(dex.args, false)  # Don't redefine
            else
                dex.args[5] = false
            end
            exprs_sigs[unwrap(dex)] = nothing
            # exprs_sigs[dex] = nothing
        end
    end
    return mod_exprs_sigs
end


text:
using PlotlyJS
include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))

using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

#D, Δx = 1.0, 0.25
D, Δx = 1.0, 0.35
Γ = meshsphere(D, Δx)
X = raviartthomas(Γ)

#Δt, Nt = 0.08, 400
Δt, Nt = 1.1, 200
T = timebasisc0d1(Δt, Nt)
U = timebasiscxd0(Δt, Nt)

V = X ⊗ T
W = X ⊗ U

width, delay, scaling = 8.0, 12.0, 1.0
gaussian = creategaussian(width, delay, scaling)

direction, polarisation = z, x
E = planewave(polarisation, direction, derive(gaussian), 1.0)
T = MWSingleLayerTDIO(1.0,-1.0,-1.0,2,0)

@hilbertspace j
@hilbertspace j′
tdefie = @discretise T[j′,j] == -1E[j′]   j∈V  j′∈W
xefie = solve(tdefie)

t1 = scatter(x=0:Δt:(Nt-1)Δt, y=xefie[1,:])

Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
ω = collect(ω0 + (0:Nt-1)*Δω)
_, i1 = findmin(abs(ω-1.0))

ω1 = ω[i1]
ue = Xefie[:,i1]/ fouriertransform(gaussian)(ω1)

fcre, geo = facecurrents(ue, X)
t2 = patch(geo, real.(norm.(fcre)))
PlotlyJS.plot(t2)


text:
@testset "abstract_gp.jl" begin
end


text:
"""
    get_map_grid_2D(par::mappingParameters)

Reconstruct the 2D grid used for mapping.
"""
function get_map_grid_2D(par::mappingParameters)

    x_grid = zeros(par.Npixels[1])
    y_grid = zeros(par.Npixels[2])

    for i = 1:par.Npixels[1]
        x_grid[i] = par.x_lim[1] + ( i - 0.5 ) * par.pixelSideLength
    end

    for i = 1:par.Npixels[2]
        y_grid[i] = par.y_lim[1] + ( i - 0.5 ) * par.pixelSideLength
    end

    return x_grid, y_grid
end

"""
    get_map_grid_3D(par::mappingParameters)

Reconstruct the 3D grid used for mapping.
"""
function get_map_grid_3D(par::mappingParameters)

    x_grid = zeros(par.Npixels[1])
    y_grid = zeros(par.Npixels[2])
    z_grid = zeros(par.Npixels[3])

    for i = 1:par.Npixels[1]
        x_grid[i] = par.x_lim[1] + ( i - 0.5 ) * par.pixelSideLength
    end

    for i = 1:par.Npixels[2]
        y_grid[i] = par.y_lim[1] + ( i - 0.5 ) * par.pixelSideLength
    end

    for i = 1:par.Npixels[3]
        z_grid[i] = par.y_lim[1] + ( i - 0.5 ) * par.pixelSideLength
    end

    return x_grid, y_grid, z_grid
end


text:
using Test
using STBParameters
using JSON
using Utils
import TBComponents: WEEKS_PER_YEAR

@testset "IT Parameter Tests" begin
    itsysdir :: IncomeTaxSys = get_default_it_system( year=2019, scotland=true)
    itsyseng :: IncomeTaxSys = get_default_it_system( year=2019, scotland=false)
    @test itsyseng.non_savings_rates[1]≈0.20
    @test itsyseng.non_savings_thresholds[2]≈150_000/WEEKS_PER_YEAR
    it = IncomeTaxSys()
    itweekly = deepcopy(it)
    @test itweekly.savings_rates == it.savings_rates
    weeklyise!( itweekly )
    annualise!( itweekly )
    @test isapprox( itweekly.non_savings_thresholds, it.non_savings_thresholds, rtol=0.00001 )
    @test itweekly.mca_minimum ≈ it.mca_minimum
    it_s = JSON.json( it )
    itj_dic = JSON.parse( it_s )
    itj = fromJSON( itj_dic )
    @test itj.non_savings_thresholds ≈ it.non_savings_thresholds
    @test itj.mca_minimum ≈ it.mca_minimum
    @test itj.company_car_charge_by_CO2_emissions ≈ it.company_car_charge_by_CO2_emissions
end # example 1


text:
struct CGrid{T <: Real, N <: Integer}
        nxc :: N
        nyc :: N
      mask :: Array{Bool,2}
         n :: N
      crop :: Diagonal{T,Array{T,1}}
      samp :: SparseMatrixCSC{T,N}
    spread :: SparseMatrixCSC{T,N}
      cent :: KronType{T,N}
end

function CGrid(;
                nxc,
                nyc,
                mask = trues(nxc,nyc))

    #check the sizes of inputs
    (size(mask) == (nxc,nyc)) || throw(DimensionMismatch("Sizes of inputs to UGrid must all be equal to nxc x nyc (i.e. $nxc x $nyc)"))

    #construct operators
    n = count(mask)
    crop = Diagonal(float(mask[:]))
    samp = crop[mask[:],:]
    spread = sparse(samp')
    cent = sparse(c(nyc)') ⊗ sparse(c(nxc)')

    #make sure boolean type rather than bitarray
    mask = convert(Array{Bool,2}, mask)

"""
CGrid(;
        nxc,
        nyc,
        mask = trues(nxc,nyc))

Construct a WAVI.jl CGrid with size (nxc,nyc)
CGrid stores operators required for the CGrid.
(Co-ordinates of CGrid stored in a Grid under xxc and yyc fields)

Keyword arguments
=================
- 'nxc': (required) Number of grid cells in x-direction in CGrid (should be same as grid.nx -1)
        Note that we store the grid size here, even though it can be easily inferred from grid, to increase transparency in velocity solve.
- 'nyc': (required) Number of grid cells in y-direction in CGrid (should be same as grid.ny - 1)
- 'mask': Mask specifying the model domain with respect to C grid
"""
return CGrid(
            nxc,
            nyc,
            mask,
            n,
            crop,
            samp,
            spread,
            cent)
end


text:
const CharString = Union{Char, AbstractString}

const NANODATE_FORMAT = dateformat"yyyy-mm-ddTHH:MM:SS.sss"

const DaysPerWeek = Int64(7)
const HoursPerDay = Int64(24)
const MinutesPerHour = Int64(60)
const SecondsPerMinute = Int64(60)
const MillisecondsPerSecond = Int64(1_000)
const MicrosecondsPerMillisecond = Int64(1_000)
const NanosecondsPerMicrosecond = Int64(1_000)

const HoursPerWeek = HoursPerDay * DaysPerWeek

const MinutesPerDay = HoursPerDay * MinutesPerHour
const MinutesPerWeek = HoursPerWeek * MinutesPerHour

const SecondsPerHour = MinutesPerHour * SecondsPerMinute
const SecondsPerDay = HoursPerDay * SecondsPerHour
const SecondsPerWeek = SecondsPerDay * DaysPerWeek

const MillisecondsPerMinute = SecondsPerMinute * MillisecondsPerSecond
const MillisecondsPerHour = MinutesPerHour * MillisecondsPerMinute
const MillisecondsPerDay = HoursPerDay * MillisecondsPerHour
const MillisecondsPerWeek = MillisecondsPerDay * DaysPerWeek

const MicrosecondsPerSecond = MillisecondsPerSecond * MicrosecondsPerMillisecond
const MicrosecondsPerMinute = SecondsPerMinute * MillisecondsPerSecond
const MicrosecondsPerHour = MinutesPerHour * MillisecondsPerMinute
const MicrosecondsPerDay = HoursPerDay * MillisecondsPerHour
const MicrosecondsPerWeek = MicrosecondsPerDay * DaysPerWeek

const NanosecondsPerMillisecond = MicrosecondsPerMillisecond * NanosecondsPerMicrosecond
const NanosecondsPerSecond = MillisecondsPerSecond * NanosecondsPerMillisecond
const NanosecondsPerMinute = SecondsPerMinute * NanosecondsPerSecond
const NanosecondsPerHour = MinutesPerHour * NanosecondsPerMinute
const NanosecondsPerDay = HoursPerDay * NanosecondsPerHour
const NanosecondsPerWeek = NanosecondsPerDay * DaysPerWeek

const Time0 = Time(0)
const Date0 = Date(0)
const DateTime0 = DateTime(0)

const Nanosecond0 = Nanosecond(0)
const Microsecond0 = Nanosecond(0)
const Millisecond0 = Millisecond(0)
const Second0 = Second(0)
const Minute0 = Minute(0)
const Hour0 = Hour(0)
const Day0 = Day(0)
const Week0 = Week(0)
const Month0 = Month(0)
const Quarter0 = Quarter(0)
const Year0 = Year(0)

const RataDie0 = Date(0, 12, 31)
const UnixDate0 = Date(1970, 1, 1)
const UnixTime0 = Time0
const UnixDateTime0 = DateTime(1970,1,1, 0,0,0,0)

const UnixRataDie = Dates.value(UnixDate0)
const Epoch2000RataDie = Dates.value(Date(2000,1,1))
const UnixMilliseconds = Dates.value(UnixDateTime0)
const UnixMicroseconds = UnixMilliseconds * MicrosecondsPerMillisecond
const UnixNanoseconds = Int128(UnixMilliseconds) * NanosecondsPerMillisecond

const Epoch2000Time = Time0
const Epoch2000Date = Date(2000,1,1)
const Epoch2000DateTime = DateTime(2000,1,1, 0,0,0, 0)
const Epoch2000RataDie = Dates.value(Epoch2000Date)
const Epoch2000Milliseconds = Dates.value(Epoch2000DateTime)
const Epoch2000Nanoseconds  = Int128(Epoch2000Milliseconds) * NanosecondsPerMillisecond


text:
for f in [:chebyshev_basis, :sin_basis, :cos_basis, :fourier_basis, :polynomial_basis, :monomial_basis]
    @eval $f(x, c) = $f(scalarize(x), c)
end
function _generateBasis!(eqs, f, x, coeffs)
    n_x = size(x, 1)
    @assert length(eqs) == size(x, 1)*length(coeffs)
    @inbounds for (i, ti) in enumerate(coeffs)
        eqs[(i-1)*n_x+1:i*n_x] .= f(x, ti)
    end
    return
end

"""
$(SIGNATURES)

Constructs an array containing a Chebyshev basis in the variables `x` with coefficients `c`.
If `c` is an `Int` returns all coefficients from 1 to `c`.
"""
function chebyshev_basis(x::Array, coefficients::AbstractVector)
    eqs = Array{Num}(undef, size(x, 1)*length(coefficients))
    f(x, t) = cos.(t .* acos.(x))
    _generateBasis!(eqs, f, x, coefficients)
    eqs
end

chebyshev_basis(x::Array, terms::Int) = chebyshev_basis(x, 1:terms)

"""
$(SIGNATURES)

Constructs an array containing a Sine basis in the variables `x` with coefficients `c`.
If `c` is an `Int` returns all coefficients from 1 to `c`.
"""
function sin_basis(x::Array, coefficients::AbstractVector)
    eqs = Array{Num}(undef, size(x, 1)*length(coefficients))
    f(x, t) = sin.(t .* x)
    _generateBasis!(eqs, f, x, coefficients)
    eqs
end

sin_basis(x::Array, terms::Int) = sin_basis(x, 1:terms)

"""
$(SIGNATURES)

Constructs an array containing a Cosine basis in the variables `x` with coefficients `c`.
If `c` is an `Int` returns all coefficients from 1 to `c`.
"""
function cos_basis(x::Array, coefficients::AbstractVector)
    eqs = Array{Num}(undef, size(x, 1)*length(coefficients))
    f(x, t) = cos.(t .* x)
    _generateBasis!(eqs, f, x, coefficients)
    eqs
end

cos_basis(x::Array, terms::Int) = cos_basis(x, 1:terms)

"""
$(SIGNATURES)

Constructs an array containing a Fourier basis in the variables `x` with (integer) coefficients `c`.
If `c` is an `Int` returns all coefficients from 1 to `c`.
"""
function fourier_basis(x::Array, coefficients::AbstractVector{Int})
    eqs = Array{Num}(undef, size(x, 1)*length(coefficients))
    f(x, t) = iseven(t) ? cos.(t .* x ./ 2) : sin.(t .* x ./2)
    _generateBasis!(eqs, f, x, coefficients)
    eqs
end

fourier_basis(x::Array, terms::Int) = fourier_basis(x, 1:terms)

"""
$(SIGNATURES)

Constructs an array containing a polynomial basis in the variables `x` up to degree `c` of the form
`[x₁, x₂, x₃, ..., x₁^1 * x₂^(c-1)]`. Mixed terms are included.
"""
function polynomial_basis(x::Array, degree::Int = 1)
    @assert degree > 0
    n_x = length(x)
    n_c = binomial(n_x+degree, degree)
    eqs = Array{Num}(undef, n_c)
    _check_degree(x) = sum(x)<=degree ? true : false
    itr = Base.Iterators.product([0:degree for i in 1:n_x]...)
    itr_ = Base.Iterators.Stateful(Base.Iterators.filter(_check_degree, itr))
    filled = false
    @inbounds for i in 1:n_c
        eqs[i] = 1
        filled = true
        for (xi, ci) in zip(x, popfirst!(itr_))
            if !iszero(ci)
                filled ? eqs[i] = xi^ci : eqs[i] *= xi^ci
                filled = false
            end
        end
    end
    eqs
end

"""
$(SIGNATURES)

Constructs an array containing monomial basis in the variables `x` up to degree `c` of the form
`[x₁, x₁^2, ... , x₁^c, x₂, x₂^2, ...]`.
"""
function monomial_basis(x::AbstractArray, degree::Int = 1)
    @assert degree > 0
    n_x = length(x)
    exponents = 1:degree
    n_e = length(exponents)
    n_c = n_x * n_e + 1
    eqs = Array{Num}(undef, n_c)
    eqs[1] = Num(1)
    idx = 0
    for i in 1:n_x, j in 1:n_e
        idx = (i-1)*n_e+j+1
        eqs[idx] = x[i]^exponents[j]
    end
    eqs
end


text:
export
    BSpline,

    NoInterp,
    Constant,
    Linear,
    Quadratic,
    Cubic

abstract Degree{N} <: Flag

immutable BSpline{D<:Degree} <: InterpolationType end
BSpline{D<:Degree}(::D) = BSpline{D}()

bsplinetype{D<:Degree}(::Type{BSpline{D}}) = D

immutable BSplineInterpolation{T,N,TCoefs<:AbstractArray,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},pad} <: AbstractInterpolation{T,N,IT,GT}
    coefs::TCoefs
end
function BSplineInterpolation{N,Tel,TWeights<:Real,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},pad}(::Type{TWeights}, A::AbstractArray{Tel,N}, ::IT, ::GT, ::Val{pad})
    isleaftype(IT) || error("The b-spline type must be a leaf type (was $IT)")
    isleaftype(typeof(A)) || warn("For performance reasons, consider using an array of a concrete type (typeof(A) == $(typeof(A)))")

    c = one(TWeights)
    for _ in 2:N
        c *= c
    end
    T = typeof(c * one(Tel))

    BSplineInterpolation{T,N,typeof(A),IT,GT,pad}(A)
end

# Utilities for working either with scalars or tuples/tuple-types
iextract{T<:BSpline}(::Type{T}, d) = T
iextract(t, d) = t.parameters[d]
padding{T,N,TCoefs,IT,GT,pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}) = pad
padding(itp::AbstractInterpolation) = padding(typeof(itp))
padextract(pad::Integer, d) = pad
padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]

lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = 1
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = size(itp, d)
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = 0.5
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = size(itp, d)+0.5

count_interp_dims{T,N,TCoefs,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType},pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}, n) = count_interp_dims(IT, n)

function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
    d <= N ? size(itp.coefs, d) - 2*padextract(pad, d) : 1
end

function interpolate{TWeights,TC,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, ::Type{TC}, A, it::IT, gt::GT)
    Apad, Pad = prefilter(TWeights, TC, A, IT, GT)
    BSplineInterpolation(TWeights, Apad, it, gt, Pad)
end
function interpolate{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, it::IT, gt::GT)
    interpolate(tweight(A), tcoef(A), A, it, gt)
end

# We can't just return a tuple-of-types due to julia #12500
tweight(A::AbstractArray) = Float64
tweight(A::AbstractArray{Float32}) = Float32
tweight(A::AbstractArray{Rational{Int}}) = Rational{Int}
tweight{T<:Integer}(A::AbstractArray{T}) = typeof(float(zero(T)))

tcoef(A::AbstractArray) = eltype(A)
tcoef(A::AbstractArray{Float32}) = Float32
tcoef(A::AbstractArray{Rational{Int}}) = Rational{Int}
tcoef{T<:Integer}(A::AbstractArray{T}) = typeof(float(zero(T)))

interpolate!{TWeights,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(::Type{TWeights}, A, it::IT, gt::GT) = BSplineInterpolation(TWeights, prefilter!(TWeights, A, IT, GT), it, gt, Val{0}())
function interpolate!{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArray, it::IT, gt::GT)
    interpolate!(tweight(A), A, it, gt)
end

include("constant.jl")
include("linear.jl")
include("quadratic.jl")
include("cubic.jl")
include("indexing.jl")
include("prefiltering.jl")
include("../filter1d.jl")


text:
using Documenter, VarianceComponentModels

ENV["DOCUMENTER_DEBUG"] = "true"
makedocs()
deploydocs(
  deps   = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
  repo   = "github.com:OpenMendel/VarianceComponentModels.jl.git",
  julia  = "0.6",
  osname = "linux"
  )


text:
# Define a "model" type, that just needs a predict function
type GenericModel
	predict # Function that makes predictions
end

type LinearModel
	predict # Funcntion that makes predictions
	w # Weight vector
end

type CompressModel
	compress # Function that compresses
	expand # Function that de-compresses
	W # weight matrix
end

# Function to compute the mode of a vector
function mode(x)
	# Returns mode of x
	# if there are multiple modes, returns the smallest
	x = sort(x[:]);

	commonVal = [];
	commonFreq = 0;
	x_prev = NaN;
	freq = 0;
	for i in 1:length(x)
		if(x[i] == x_prev)
			freq += 1;
		else
			freq = 1;
		end
		if(freq > commonFreq)
			commonFreq = freq;
			commonVal = x[i];
		end
		x_prev = x[i];
	end
	return commonVal
end

# Return element-wise log, but set log(0)=0
function log0(x)
	y = copy(x)
	y[y.==0] = 1
	return log.(y)
end

# Return squared Euclidean distance all pairs of rows in X1 and X2
function distancesSquared(X1,X2)
	(n,d) = size(X1)
	(t,d2) = size(X2)
	assert(d==d2)
	return X1.^2*ones(d,t) + ones(n,d)*(X2').^2 - 2X1*X2'
end

# Subtract mean of each column and divide by standard deviation
# (or call it with mu and sigma to use these specific mean/std)
function standardizeCols(X;mu=[],sigma=[])
	(n,d) = size(X)

	if isempty(mu)
		mu_j = mean(X,1)
	else
		mu_j = mu
	end

	Xstd = zeros(n,d)
	for j in 1:d
		Xstd[:,j] = X[:,j] - mu_j[j]
	end

	if isempty(sigma)
		sigma_j = std(Xstd,1)
	else
		sigma_j = sigma
	end

	for j in 1:d
		Xstd[:,j] /= sigma_j[j]
	end

	if isempty(mu) & isempty(sigma)
		return (Xstd,mu_j,sigma_j)
	else
		return Xstd
	end
end

### A function to compute the gradient numerically
function numGrad(func,x)
	n = length(x);
	delta = 2*sqrt(1e-12)*(1+norm(x));
	g = zeros(n);
	e_i = zeros(n)
	for i = 1:n
		e_i[i] = 1;
		(fxp,) = func(x + delta*e_i)
		(fxm,) = func(x - delta*e_i)
		g[i] = (fxp - fxm)/2delta;
		e_i[i] = 0
	end
	return g
end

### Check if number is a real-finite number
function isfinitereal(x)
	return (imag(x) == 0) & (!isnan(x)) & (!isinf(x))
end

function dijkstra(D,s,t)
	# Finds shortest path from 's' to 't' with non-negative edge weights D
  # - set the edge weight D[i,j] to Inf if there is no edge between nodes 'i' and 'j' in the graph
	# (not really a fast implementation)

	n = size(D,1)

	distances = fill(Inf,n)
	distances[s] = 0

	visited = fill(false,n)

	while true

		# Choose next node to visit
		minDist = Inf
		i = []
		for j in 1:n
			if !visited[j] & (distances[j] < minDist)
				i = j
				minDist = distances[j]
			end
		end

		if isinf(minDist)
			# We have nowhere left to visit, there is no path
			break
		end

		# "Visit" the node: update neighbour distances
		for j in 1:n
			if distances[i] + D[i,j] < distances[j]
				distances[j] = distances[i] + D[i,j]
			end
		end
		visited[i] = true

		# If we just visited the solution
		if i == t
			break
		end
	end
	return distances[t]
end


text:
using Documenter, DocumenterTools, DocumenterLaTeX
using Test

const ROOT = joinpath(@__DIR__, "..")

# Documenter package docs
doc = makedocs(
    debug = true,
    root = ROOT,
    build = "pdf/build",
    modules = [Documenter, DocumenterTools],
    clean = false,
    format = LaTeX(platform = "docker"),
    sitename = "Documenter.jl",
    authors = "Michael Hatherly, Morten Piibeleht, and contributors.",
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/examples.md",
            "man/syntax.md",
            "man/doctests.md",
            "man/latex.md",
            hide("man/hosting.md", [
                "man/hosting/walkthrough.md"
            ]),
            "man/other-formats.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            hide("Internals" => "lib/internals.md", Any[
                "lib/internals/anchors.md",
                "lib/internals/builder.md",
                "lib/internals/cross-references.md",
                "lib/internals/docchecks.md",
                "lib/internals/docsystem.md",
                "lib/internals/doctests.md",
                "lib/internals/documenter.md",
                "lib/internals/documentertools.md",
                "lib/internals/documents.md",
                "lib/internals/dom.md",
                "lib/internals/expanders.md",
                "lib/internals/mdflatten.md",
                "lib/internals/selectors.md",
                "lib/internals/textdiff.md",
                "lib/internals/utilities.md",
                "lib/internals/writers.md",
            ])
        ],
        "contributing.md",
    ]
);

# hack to only deploy the actual pdf-file
mkpath(joinpath(ROOT, "pdf", "build", "pdfdir"))
let files = readdir(joinpath(ROOT, "pdf", "build"))
    for f in files
        if startswith(f, "Documenter.jl") && endswith(f, ".pdf")
            mv(joinpath(ROOT, "pdf", "build", f),
               joinpath(ROOT, "pdf", "build", "pdfdir", f))
        end
    end
end

deploydocs(
    repo = "github.com/JuliaDocs/Documenter.jl.git",
    root = ROOT,
    target = "pdf/build/pdfdir",
    branch = "gh-pages-pdf",
    forcepush = true,
)


text:
# getSHbasis.jl
using SphericalHarmonics: computeYlm
using CoordinateTransformations: SphericalFromCartesian
using MIRTjim: jim

# export numSH getSHbasis

"""
    c2sph(r, l, m)

Evaluate the `(l,m)` spherical harmonic at one spatial location `(x, y, z)`.
"""
function c2sph(x, y, z, l, m)
	a = SphericalFromCartesian()([x,y,z]) # todo: SVector?
	(rho, ϕ, θ) = (a.r, a.θ, π/2 - a.ϕ)   # (radius, azimuth, colatitude)
	Y = computeYlm(θ, ϕ; lmax=l)
	rho^l * Y[(l,m)]
end

"""
    n = numSH(L::Int)
Number of spherical harmonics for given `L`
i.e., `sum(2*(0:L) .+ 1)`, which is 9 for `L=2`.
Order of 0th-2nd order terms (for `L=2`):
* 1 cf (center frequency, Hz)
* 2 z
* 3 x
* 4 y
* 5 z2
* 6 zx
* 7 zy
* 8 x2y2
* 9 xy
"""
function numSH(L::Int)
	sum(2*(0:L) .+ 1)
end

"""
    getSHbasis!(h, x, y, z; L=2)

Store in `h` the spherical harmonic basis up to order `L` (default 2),
evaluated at spatial location `(x, y, z)`.
Vector `h` must have length `numSH(L)`.
"""
function getSHbasis!(h::AbstractVector{<:Real}, x::Real, y::Real, z::Real; L::Int = 2)

	ic = 1
	for l = 0:L, m = 0:l
		f = c2sph(x, y, z, l, m)
		h[ic] = real(f)
		ic += 1
		if m != 0
			h[ic] = imag(f)
			ic += 1
 		end
 	end

	h[1] = 1.0 # basis for center frequency offset

	return h
end

"""
    H = getSHbasis(x, y, z; L::Int=2)

Get spherical harmonic basis up to order `L` (default 2)
evaluated at spatial location vectors `x, y, z`
See `numSH` for ordering.
Output size is `(length(x), length(y), length(z), numSH(L))`
"""
function getSHbasis(
	x::AbstractVector{<:Real},
	y::AbstractVector{<:Real},
	z::AbstractVector{<:Real};
	L::Int=2
)

	T = promote_type(eltype.([x, y, z])..., Float32) # at least Float32
	H = zeros(T, length(x), length(y), length(z), numSH(L))
	for ix=1:length(x), iy=1:length(y), iz=1:length(z)
		getSHbasis!((@view H[ix,iy,iz,:]), x[ix], y[iy], z[iz]; L)
	end
	return H
end

# convenience method for scalar inputs
getSHbasis(x::Real, y::Real, z::Real; L::Int=2) = getSHbasis([x], [y], [z]; L)

"""
    H = getSHbasis("test"; L::Int = 2)

Test function
"""
function getSHbasis(str::String; L::Int = 2,
	nx::Int=22, ny::Int=20, nz::Int=18, fov::NTuple{3,Real}=(20,20,20),
)
	rx = LinRange(-1,1,nx)*fov[1]/2
	ry = LinRange(-1,1,ny)*fov[2]/2
	rz = LinRange(-1,1,nz)*fov[3]/2
	@time H = getSHbasis(rx, ry, rz; L)
	jim(H; ncol=numSH(L), color=:jet, gui=true, line3plot=false) # compare with >> evalspharm("test")
	return H
end

# getSHbasis("test"; L = 2);


text:
module TestOpenCL
using Test
using OpenCL
using Base.GC

@testset "layout" begin
    x = ((10f0, 1f0, 2f0), (10f0, 1f0, 2f0), (10f0, 1f0, 2f0))
    clx = cl.replace_different_layout(x)

    @test clx == ((10f0, 1f0, 2f0, 0f0), (10f0, 1f0, 2f0, 0f0), (10f0, 1f0, 2f0, 0f0))
    x = (nothing, nothing, nothing)
    clx = cl.replace_different_layout(x)
    @test clx == 0 # TODO should it be like this?
end

function create_test_buffer()
    ctx = cl.create_some_context()
    queue = cl.CmdQueue(ctx)
    testarray = zeros(Float32, 1000)
    buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=testarray)
    return (queue, buf, testarray, ctx)
end

include("test_platform.jl")
include("test_context.jl")
include("test_device.jl")
include("test_cmdqueue.jl")
include("test_minver.jl")
include("test_event.jl")
include("test_program.jl")
include("test_kernel.jl")
include("test_behaviour.jl")
include("test_memory.jl")
include("test_buffer.jl")
include("test_array.jl")

@testset "context jl reference counting" begin
    Base.GC.gc()
    @test isempty(cl._ctx_reference_count)
end

end # module


text:
using Test

@testset "All tests" begin
    include("earth_models_1d.jl")
    include("inbuilt_models.jl")
    include("basic_properties.jl")
end


text:
#!/usr/bin/env julia

using DatagenCopulaBased
import DatagenCopulaBased: ρ2θ
using Distributions
using Cumulants
using PyPlot
using SymmetricTensors
using StatsBase
using PyCall
using JLD2
using FileIO
using Distributed
using HCubature

#using QuadGK
#using Cubature
#using NLsolve
#using LsqFit

mpl = pyimport("matplotlib")
mpl.rc("text", usetex=true)

addprocs(3)
@everywhere using DatagenCopulaBased
@everywhere using Cumulants
nprocs()

function ccl(u, θ::Float64)
  q = length(u)
  (1-q+sum(u.^(-θ)))^(-q-1/θ)*prod([u[j]^(-θ-1)*((j-1)*θ+1) for j in 1:q])
end

function munif(i::Vector{Int}, θ::Union{Int, Float64}, g::Function)
  n = length(i)
  f(x) = mapreduce(j -> x[j], *, i)*g(x, θ)
  hcubature(f, [fill(0, n)...],[fill(1, n)...])[1]
end

function c3(θ::Union{Int, Float64}, g::Function)
 c3ofd = munif([1,2,3], θ, g)-3*munif([1,2], θ, g)*munif([1], θ, g)+2*munif([1], θ, g)^3
 c3mix = munif([1,1,2], θ, g)-munif([1,1], θ, g)*munif([1], θ, g)-2*munif([1,2], θ, g)*munif([1], θ, g)+2*munif([1], θ, g)^3
 c3diag = munif([1,1,1], θ, g)-3*munif([1,1], θ, g)*munif([1], θ, g)+2*munif([1], θ, g)^3
 [c3diag, c3mix, c3ofd]
end

function mmarg(i::Vector{Int}, θ::Union{Int, Float64}, g::Function, d)
  n = length(i)
  f(x) = mapreduce(j -> quantile.(d, x[j]), *, i)*g(x, θ)
  hcubature(f, [fill(0.00001, n)...],[fill(0.999999, n)...])[1]
end

c3symmarg(θ::Union{Int, Float64}, g::Function, d) =[mmarg([1,1,1], θ, g, d), mmarg([1,1,2], θ, g, d), mmarg([1,2,3], θ, g, d)]

function c4symmarg(θ::Union{Int, Float64}, g::Function, d)
  a = mmarg([1,1,1,1], θ, g, d) - 3*mmarg([1,1], θ, g, d)^2
  println("a")
  b = mmarg([1,1,1,2], θ, g, d) - 3*mmarg([1,1], θ, g, d)*mmarg([1,2], θ, g, d)
  println("b")
  c = mmarg([1,1,2,2], θ, g, d) - mmarg([1,1], θ, g, d)^2 - 2*mmarg([1,2], θ, g, d)^2
  println("c")
  e = mmarg([1,1,2,3], θ, g, d) - mmarg([1,1], θ, g, d)*mmarg([1,2], θ, g, d) - 2*mmarg([1,2], θ, g, d)^2
  println("d")
  f = mmarg([1,2,3,4], θ, g, d) - 3*mmarg([1,2], θ, g, d)^2
  [a,b,c,e,f]
end

function c3empirical(cop, d, t::Int = 1_000_000)
  u = quantile.(d, simulate_copula(t, cop))
  c = cumulants(u, 3)[3]
  [c[1,1,1], c[1,1,2], c[1,2,3]]
end

function c4empirical(cop, d, t::Int = 5000000)
  u = quantile.(d, simulate_copula(t, cop))
  c = cumulants(u, 4)[4]
  [c[1,1,1,1], c[1,1,1,2], c[1,1,2,2], c[1,1,2,3], c[1,2,3,4]]
end

function stats3theoreunifcl()
  θ = map(i -> ρ2θ(i, "clayton"), 0.05:0.05:.95)
  n = length(θ)
  t = zeros(n, 3)
  for i in 1:n
    t[i,:] = c3(θ[i], ccl)
  end
  t
end

function statstheoreticalclayton(m::Int = 3, d = Normal(0,1), ρ)
  θ = map(i -> ρ2θ(i, "clayton"), ρ)
  n = length(θ)
  t = (m == 3) ? zeros(n, 3) : zeros(n, 5)
  for i in 1:n
    t[i,:] = (m == 3) ? c3symmarg(θ[i], ccl, d) : c4symmarg(θ[i], ccl, d)
    println(i)
  end
  t
end

#save("c3clgmarg.jld2", "emp", e, "theoret", t, "rho", ρ)

function empiricalcums(copula, m::Int, ρ::Vector{Float64})
  println(ρ)
  n = length(ρ)
  e = (m == 3) ? zeros(n, 3) : zeros(n, 5)
  for i in 1:n
    println(i)
    e[i,:] = (m == 3) ? c3empirical(copula(3, ρ[i], "Spearman"), Normal(0,1)) : c4empirical(copula(4, ρ[i], "Spearman"), Normal(0,1))
  end
  e
end

function plotc3empth(e,t, ρ, ρ1)
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  fx = matplotlib.ticker.ScalarFormatter()
  fx.set_powerlimits((-1, 4))
  ax.yaxis.set_major_formatter(fx)
  plot(ρ, e[:,3], "o", color = "black", label = "\$ \\mathbf{i} = (1,2,3) \$", markersize = 1.)
  plot(ρ1, t[:,3], color = "black", linewidth = 0.5)
  plot(ρ, e[:,2], "o", color = "blue", label = "\$ \\mathbf{i} = (1,1,2) \$", markersize = 1.)
  plot(ρ1, t[:,2], color = "blue", linewidth = 0.5)
  plot(ρ, e[:,1], "o", color = "red", label = "\$ \\mathbf{i} = (1,1,1) \$", markersize = 1.)
  plot(ρ1, t[:,1], color = "red", linewidth = 0.5)
  PyPlot.ylabel("cumulant element", labelpad = -1.0)
  PyPlot.xlabel("Spearman \$ \\rho \$ between marginals", labelpad = -1.0)
  ax.legend(fontsize = 6, loc = 3, ncol = 1)
  subplots_adjust(left = 0.12, bottom = 0.12)
  fig.savefig("pics/c3clgmarg.pdf")
end

function plotc3emp(e, ρ, cop::String)
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  fx = matplotlib.ticker.ScalarFormatter()
  fx.set_powerlimits((-1, 4))
  ax.yaxis.set_major_formatter(fx)
  plot(ρ, e[:,3], "o", color = "black", label = "\$ \\mathbf{i} = (1, 2, 3) \$", markersize = 1.)
  plot(ρ, e[:,2], "o", color = "blue", label = "\$ \\mathbf{i} = (1, 1, 2) \$", markersize = 1.)
  plot(ρ, e[:,1], "o", color = "red", label = "\$ \\mathbf{i} = (1, 1, 1) \$", markersize = 1.)
  PyPlot.ylabel("cumulant element", labelpad = -1.0)
  PyPlot.xlabel("Spearman \$ \\rho \$ between marginals", labelpad = -1.0)
  #ax.legend(fontsize = 7, loc = 3, ncol = 1)
  subplots_adjust(left = 0.16, bottom = 0.12)
  fig.savefig("pics/c3"*cop*"gmarg.pdf")
end

function plotc4emp(e, ρ, cop::String)
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  fx = matplotlib.ticker.ScalarFormatter()
  fx.set_powerlimits((-1, 4))
  ax.yaxis.set_major_formatter(fx)
  plot(ρ, e[:,5], "o", color = "black", label = "\$ i_1, i_2, i_3, i_4 \$", markersize = 1.)
  plot(ρ, e[:,4], "o", color = "brown", label = "\$ i_1, i_1, i_2, i_3 \$", markersize = 1.)
  plot(ρ, e[:,3], "o", color = "gray", label = "\$ i_1, i_1, i_2, i_2 \$", markersize = 1.)
  plot(ρ, e[:,2], "o", color = "blue", label = "\$ i_1, i_1, i_1, i_2 \$", markersize = 1.)
  plot(ρ, e[:,1], "o", color = "red", label = "\$ i_1, i_1, i_1, i_1 \$", markersize = 1.)
  PyPlot.ylabel("cumulant element", labelpad = -1.0)
  PyPlot.xlabel("Spearman \$ \\rho \$", labelpad = -1.0)
  #ax.legend(fontsize = 4.5, loc = 3, ncol = 1)
  subplots_adjust(left = 0.16, bottom = 0.12)
  fig.savefig("pics/c4"*cop*"gmarg.pdf")
end

function plotc4empth(e,t, ρ, ρ1)
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  fx = matplotlib.ticker.ScalarFormatter()
  fx.set_powerlimits((-1, 4))
  ax.yaxis.set_major_formatter(fx)
  plot(ρ1, t[:,5], color = "black", linewidth = 0.5)
  plot(ρ, e[:,5], "o", color = "black", label = "\$ i_1, i_2, i_3, i_4 \$", markersize = 1.)
  plot(ρ1, t[:,4], color = "brown", linewidth = 0.5)
  plot(ρ, e[:,4], "o", color = "brown", label = "\$ i_1, i_1, i_2, i_3 \$", markersize = 1.)
  plot(ρ, e[:,3], "o", color = "gray", label = "\$ i_1, i_1, i_2, i_2 \$", markersize = 1.)
  plot(ρ1, t[:,3], color = "gray", linewidth = 0.5)
  plot(ρ1, t[:,2], color = "blue", linewidth = 0.5)
  plot(ρ, e[:,2], "o", color = "blue", label = "\$ i_1, i_1, i_1, i_2 \$", markersize = 1.)
  plot(ρ, e[:,1], "o", color = "red", label = "\$ i_1, i_1, i_1, i_1 \$", markersize = 1.)
  plot(ρ1, t[:,1], color = "red", linewidth = 0.5)
  PyPlot.ylabel("cumulant element", labelpad = -1.0)
  PyPlot.xlabel("Spearman \$ \\rho \$", labelpad = -1.0)
  ax.legend(fontsize = 4.5, loc = 3, ncol = 1)
  subplots_adjust(left = 0.12, bottom = 0.12)
  fig.savefig("pics/c4clgmarg.pdf")
end

if false
  ρ = collect(0.05:0.05:.95)
  e4 = empiricalcums(Clayton_cop, 4, ρ)
  save("pics/c4clgmarg.jld2", "emp", e4,"rho", ρ)
  t = statstheoreticalclayton(4, Normal(0,1), ρ);
  save("pics/teorc4clgmarg.jld2", "theor", t,"rho", ρ)
end
if false
  ρ = collect(0.05:0.05:.95)
  e3 = empiricalcums(Clayton_cop, 3, ρ)
  save("pics/c3clgmarg.jld2", "emp", e3,"rho", ρ)
  t = statstheoreticalclayton(3, Normal(0,1), ρ);
  save("pics/teorc3clgmarg.jld2", "theor", t,"rho", ρ)
end
if false
  ρ = [i for i in 0.02:0.02:.85]
  e = empiricalcums(Frank_cop, 3, ρ)
  save("pics/c3frgmarg.jld2", "emp", e,"rho", ρ)
end
if false
  ρ = [i for i in 0.02:0.02:.98]
  e = empiricalcums(Gumbel_cop, 3, ρ)
  save("pics/c3gugmarg.jld2", "emp", e,"rho", ρ)
end
if false
  ρ = vcat([i for i in 0.02:0.02:.49], [0.499])
  e = empiricalcums(AMH_cop, 3, ρ)
  save("pics/c3amhgmarg.jld2", "emp", e,"rho", ρ)
end

emp = load("pics/c3clgmarg.jld2")
t = load("pics/teorc3clgmarg.jld2")
plotc3empth(emp["emp"], t["theor"], emp["rho"], t["rho"])

c3els = load("pics/c3gugmarg.jld2")
plotc3emp(c3els["emp"], c3els["rho"], "gu")

c3els = load("pics/c3frgmarg.jld2")
plotc3emp(c3els["emp"], c3els["rho"], "fr")

c3els = load("pics/c3amhgmarg.jld2")
plotc3emp(c3els["emp"], c3els["rho"], "amh")

emp = load("pics/c4clgmarg.jld2")
t = load("pics/teorc4clgmarg.jld2")
plotc4empth(emp["emp"], t["theor"], emp["rho"], t["rho"])


text:
"""
    Maybe{T} = Union{Nothing, T}

A nice name for "nullable" type.
"""
const Maybe{T} = Union{Nothing,T}

"""
    _default(d::T, x::Maybe{T})::T where {T}

Fold the `Maybe{T}` down to `T` by defaulting.
"""
function _default(d::T, x::Maybe{T})::T where {T}
    isnothing(x) ? d : x
end

"""
    _maybemap(f, x::Maybe)::Maybe

Apply a function to `x` only if it is not `nothing`.
"""
function _maybemap(f, x::Maybe)::Maybe
    isnothing(x) ? nothing : f(x)
end


text:
"""
    genericTaskCode
This is a template module for keeping records of day to day experiments.

# This scripts contains the following tests:

# `\$DIR` this file supposed to run in

# task result summary

# see also

# tags

"""
module genericTaskCode

"""
    run
The script main content.
"""
function run()
  println("Blank script.")
  nothing
end

end


text:
using Test
using StaticArrays
using Proj4
import PROJ_jll

function read_cmd(cmd)
    bytes = read(cmd)
    str = String(bytes)
    # on Windows there are CFLR line endings
    lf_only = replace(str, "\r" => "")
    no_tabs = replace(lf_only, "\t" => " ")
    return String(strip(no_tabs))
end

# https://proj.org/apps/cs2cs.html#using-epsg-crs-codes
@test PROJ_jll.cs2cs() do cs2cs
    read_cmd(pipeline(IOBuffer("45N 2E"), `$cs2cs EPSG:4326 EPSG:32631`))
end == "421184.70 4983436.77 0.00"

# https://proj.org/usage/transformation.html?#grid-based-datum-adjustments
@test PROJ_jll.cs2cs() do cs2cs
    withenv("PROJ_NETWORK" => "ON") do
        # instead of ntv1_can.dat we can also specify ca_nrc_ntv1_can.tif
        read_cmd(pipeline(IOBuffer("-111 50"),
        `$cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat
            +to +proj=latlong +ellps=GRS80 +datum=NAD83`))
    end
end == "111d0'3.006\"W 50d0'0.125\"N 0.000"

@test PROJ_jll.cs2cs() do cs2cs
    withenv("PROJ_NETWORK" => "ON") do
        read_cmd(pipeline(IOBuffer("-111 50"),
        `$cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ca_nrc_ntv1_can.tif
            +to +proj=latlong +ellps=GRS80 +datum=NAD83`))
    end
end == "111d0'3.006\"W 50d0'0.125\"N 0.000"

function xyzt_transform_cli(point::AbstractVector; network::Bool=false)
    # convert the vector to a string like "-33 151 5 2020"
    input = String(strip(string(point'), ('[', ']')))
    xyzt_transform_cli(input; network)
end

function xyzt_transform_cli(point::String; network::Bool=false)
    proj_network = network ? "ON" : nothing
    PROJ_jll.cs2cs() do cs2cs
        withenv("PROJ_NETWORK" => proj_network) do
            read_cmd(pipeline(IOBuffer(point),
                `$cs2cs -d 6 EPSG:4326+5773 EPSG:7856+5711`))
        end
    end
end

function xyzt_transform(point::AbstractVector; network::Bool=false)
    proj_network = network ? "ON" : nothing
    trans = Proj4.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true)
    trans(point)
end

# http://osgeo-org.1560.x6.nabble.com/PROJ-Different-results-from-cs2cs-vs-C-API-code-td5446396.html

# When you use cs2cs and do not specify it (as the 4th value in the coordinate),
# the time is sent to HUGE_VAL, which makes it not being used at all
# (so here as the reference epoch for the time-dependent Helmert transformation is 2020.0,
# it will be as if you specified 2020.0).

xy1 = "313152.777216 6346936.495714"
xy2 = "313152.777214 6346936.495810"
@test xyzt_transform_cli(SA[-33, 151, 5, 2020]; network=true) == "$xy1 5.280678 2020"
@test xyzt_transform_cli("-33 151 5 2020"; network=true) == "$xy1 5.280678 2020"
# no z correction with no network, as expected
@test xyzt_transform_cli(SA[-33, 151, 5, 2020]; network=false) == "$xy2 5.000000 2020"
# year 0, different xyz, as expected
@test xyzt_transform_cli(SA[-33, 151, 5, 0]; network=true) == "313188.897461 6347047.288591 4.940935 0"
# no time input is like passing 2020, as expected
@test xyzt_transform_cli(SA[-33, 151, 5]; network=true) == "$xy1 5.280678"

# need to switch axis order here
# interestingly this does not respond to the PROJ_NETWORK environment variable,
# but does respond to enabling it in the context
Proj4.proj_context_is_network_enabled()
@test Proj4.proj_context_set_enable_network(1) == 1
# this is like passing floatmax() to the cli, but if we do it correctly we expect
# it to be like the version that passes only xyz to the cli
@test xyzt_transform(SA_F64[151, -33, 5]) == SA[313152.7772137531, 6.346936495809965e6, 5.280647277836724]
# this is expected, like the above
@test xyzt_transform(SA_F64[151, -33, 5, 2020]) == SA[313152.77721557155, 6.34693649571435e6, 5.28067830334755, 2020.0]

@test xyzt_transform(SA_F64[151, -33, 5, floatmax() * 2]) == SA[313152.77721557155, 6.34693649571435e6, 5.28067830334755, Inf]

@test PROJ_jll.cs2cs() do cs2cs
    read_cmd(pipeline(IOBuffer("-33 151 5 0"), `$cs2cs -d 8 EPSG:4326+5773 EPSG:7856+5711`))
end == "313152.77721375 6346936.49580996 5.00000000 0"

@test PROJ_jll.projinfo() do projinfo
    read_cmd(`$projinfo  -o PROJ EPSG:25832`)
end == """PROJ.4 string:
+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs
"""

PROJ_jll.projinfo() do projinfo
    n_candidate_ops = "Candidate operations found: 4"
    grid_not_found = "Grid us_nga_egm96_15.tif needed but not found on the system." *
        " Can be obtained at https://cdn.proj.org/us_nga_egm96_15.tif"
    info = read_cmd(`$projinfo -s EPSG:4326+5773 -t EPSG:7856+5711`)
    @test occursin(n_candidate_ops, info)
    @test occursin(grid_not_found, info)
end


text:
"""
Make a cartesian product equispaced grid on `[-1,1]^N`.
"""
function equispaced_grid(::Val{N}, n::Int) where {N}
    grid1 = PeriodicEquispacedGrid{T}(n)
    grid = ProductGrid( ntuple(t->grid1, Val{N}())... )
end

make_centers(grid::AbstractGrid) = collect(grid)

restrict_grid(grid, domain) = MaskedGrid(grid, domain)

function rbfapprox(rbf::RBF, centers, points, f; verbose = false, method = :svd)
    basis = RBFDictionary(centers, rbf)
    A = evaluation_matrix(basis, points)
    B = f.(points)
    if method == :LAPACK
        coef, rank = Base.LinAlg.LAPACK.gelsy!(A,B)
        if verbose
            println("Effective rank: ", rank)
        end
    elseif method == :svd
        threshold = 1e-14
        usv = LAPACK.gesdd!('S',A)
        s = usv[2]
        maxindex = findlast(s.>threshold)
        if verbose
            println("Effective svd-rank: ", maxindex)
        end
        sreg = 1 ./ s[1:maxindex]
        ureg = usv[1][:,1:maxindex]
        vreg = usv[3][1:maxindex,:]
        coef = vreg' * diagm(sreg) * (ureg'*B)
    else
        error("Please supply a valid method for rbfapprox")
    end

    Expansion(basis, coef)
end


text:
using SparseArrays
using LinearAlgebra
using JuMP

include("proof_helper.jl")
include("../src/hypergraph-affinity-functions.jl")
k = 3

r = round(Int64,(k+1)/2)
alph = 0.56
b1 = arbitrary_baselines(k,alph)
b2 = arbitrary_baselines(k,1-alph)

n = 100
n1 = round(Int64,floor(n*alph))
n2 = n - n1
b1, b2 = n1_n2_baselines(k,n1,n2)

Dk, Db, R, B = Bmat(b1,b2)
Xopt, Yopt, gam = PrimalDualXY(gurobi_env,B,k)
Ydual = OptimalDual(b1,b2)
Ybal, delta =  BalancedDual(k)
@assert(norm(Ydual-Yopt) < 1e-15)

## Now we start to check how things hold

V = R'*Db*Yopt
y = y_double(b1,b2)
t = k-1
c = 2

## I know how to recover V
for c = 1:2
    for t = k:-1:r
        r2t = t*sum(b2[i]*y[2,i] for i = r:k) + (k-t)*sum(b1[i]*y[1,i] for i = r:k)
        V2t = V[tt(2,t)]

        @show norm(r2t- V2t)
    end
end

## Now let's manipulate terms
t = k-1
c = 2
yvec[2,t]*t
y2k = Yopt[1]
V2t = V[tt(2,t)]

# ID (1)
s1 = sum(b2[i] for i = r:k) + sum(b1[i] for i = r:k) - b2[r]

C1 = y2k*t*sum((2-k/i)*b2[i] for i = r:k) + y2k*t*sum((2-k/i)*b1[i] for i = r:k)
C2 = y2k*t*sum((2-k/i)*(b2[i]+b1[i]) for i = r:k)

AA = 2*k*t/delta*sum(b2[i]/i for i = r:k)
BB = -2*t/delta*sum(b2[i] for i = r:k)
CC = t*y2k*sum((2-k/i)*b2[i] for i = r:k)
DD = 2*(k-t)/delta*sum(b1[i] for i = r:k)
FF = -k*y2k*sum((2-k/i)*b1[i] for i = r:k)
GG = t*y2k*sum((2-k/i)*b1[i] for i = r:k)
ans1 = AA+BB+CC+DD+FF+GG

## Checking baselines

i = k-1

k/i*b1[i]

(1-alph)^(k-i)*alph^(i-1)*binomial(k,i)

## Four piece
y = y2k
bet2 = sum(b2[i] for i = r:k)
bet1 = sum(b1[i] for i = r:k)
bet2i = sum(b2[i]/i for i = r:k)
bet1i = sum(b1[i]/i for i = r:k)
ans2 = bet2*( -2*t/delta + 2*y*t) + bet2i*(2*k*t/delta - t*y*k) + bet1*(2*(k-t)/delta + (t-k)*y*2) + bet1i*((k-t)*k*y)

## New approach GOOD (1) in notes

ls = (k-t)*sum(b1[i]*(2/delta - (2-k/i)*y) for i = r:k)
rs = (2/delta*(k - t) + (2*t-k)*y) - t*y
r2 = 2/delta*(k-t) + t*y - k*y
r3 = (2/delta - y)*(k-t)

## simplifies to...

ls = sum(b1[i]*(2/delta - (2-k/i)*y) for i = r:k)
rs = (2/delta - y)
r2 = 2/delta*( (1- sum(b2[i] for i = r:k))/( 1 - sum(b2[i]*(2-k/i) for i = r:k)))
l2 = sum(b1[i] for i = r:k)*(2/delta - 2*y) + sum(b1[i]*(k/i)*y for i = r:k)

yguess = 2/delta*(1- sum(b1[i] for i = r:k))/(1-sum(b1[i]*(2-k/i) for i = r:k))

norm(y-yguess)

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

# # Multiple neural network layers with `Flux.jl`

#-

# In a previous notebook, we saw that one layer of neurons wasn't enough to distinguish between three types of fruit (apples, bananas *and* grapes), since the data is quite complex. To solve this problem, we need to use more layers, so heading into the territory of **deep learning**!
#
# By adding another layer between the inputs and the output neurons, a so-called "hidden layer", we will get our first serious **neural network**, looking something like this:

include(datapath("scripts/draw_neural_net.jl"))
draw_network([2, 4, 3])

# We will continue to use two input data and try to classify into three types, so we will have three output neurons. We have chosen to add a single "hidden layer" in between, and have arbitrarily chosen to put 4 neurons there.
#
# Much of the *art* of deep learning is choosing a suitable structure for the neural network that will allow the model to be sufficiently complex to model the data well, but sufficiently simple to allow the parameters to be learned in a reasonable length of time.

#-

# ## Read in and process data
#
# As before, let's load some pre-processed data using code we've seen in the previous notebook.

## using Pkg; Pkg.add("Flux")
using Flux
using Flux: onehot

#-

## using Pkg; Pkg.add("CSV")
using CSV, DataFrames

apples_1 = DataFrame(CSV.File(datapath("data/Apple_Golden_1.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples_2 = DataFrame(CSV.File(datapath("data/Apple_Golden_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples_3 = DataFrame(CSV.File(datapath("data/Apple_Golden_3.dat"), delim='\t', allowmissing=:none, normalizenames=true))
bananas = DataFrame(CSV.File(datapath("data/Banana.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes_1 = DataFrame(CSV.File(datapath("data/Grape_White.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes_2 = DataFrame(CSV.File(datapath("data/Grape_White_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))

apples = vcat(apples_1, apples_2, apples_3)
grapes = vcat(grapes_1, grapes_2);

#-

col1 = :red
col2 = :blue

x_apples  = [ [apples_1[i, col1], apples_1[i, col2]] for i in 1:size(apples_1)[1] ]
append!(x_apples, [ [apples_2[i, col1], apples_2[i, col2]] for i in 1:size(apples_2)[1] ])
append!(x_apples, [ [apples_3[i, col1], apples_3[i, col2]] for i in 1:size(apples_3)[1] ])

x_bananas = [ [bananas[i, col1], bananas[i, col2]] for i in 1:size(bananas)[1] ]

x_grapes = [ [grapes_1[i, col1], grapes_1[i, col2]] for i in 1:size(grapes_1)[1] ]
append!(x_grapes, [ [grapes_2[i, col1], grapes_2[i, col2]] for i in 1:size(grapes_2)[1] ])

xs = vcat(x_apples, x_bananas, x_grapes);

# We now we wish to classify the three types of fruit, so we again use one-hot vectors to represent the desired outputs $y^{(i)}$:

labels = [ones(length(x_apples)); 2*ones(length(x_bananas)); 3*ones(length(x_grapes))];

ys = [onehot(label, 1:3) for label in labels];  # onehotbatch(labels, 1:3)

# The input data is in `xs` and the one-hot vectors are in `ys`.

#-

# ## Multiple layers in Flux

#-

# Let's tell Flux what structure we want the network to have. We first specify the number of neurons in each layer, and then construct each layer as a `Dense` layer:

inputs = 2
hidden = 4
outputs = 3

layer1 = Dense(inputs, hidden, σ)
layer2 = Dense(hidden, outputs, σ)

# To stitch together multiple layers to make a multi-layer network, we use Flux's `Chain` function:

const model = Chain(layer1, layer2)

# #### Exercise 1
#
# What is the internal structure and sub-structure of this `model` object?

#-

# #### Solution

model.layers

# So `model` understands that it consists of two layers, and:

model.layers[1].W

#-

model.layers[1].b

# So both the layers have their own `W` and `b`, just as you might expect -- `Flux` wraps all this up in a way that is structured but easy to use.
#
# In particular, using `params` returns *all* the parameters hidden inside the `model` object, that is the pairs $(W, b)$ from both layers:

params(model)

# ## Training the model

#-

# We have now set up a model and we have some training data.
# How do we train the model on the data?
#
# The amazing thing is that the rest of the code in `Flux` is **exactly the same as before**. This is possible thanks to the design of Julia itself, and of the `Flux` package.

#-

# #### Exercise 2
#
# Train the model as before, now using the popular `ADAM` optimizer. You may need to train the network for longer than before, since we have many more parameters.

#-

# #### Solution

loss(x, y) = Flux.mse(model(x), y)

#-

data = zip(xs, ys);

#-

opt = ADAM(params(model), 0.02)

#-

@time for i in 1:100
    Flux.train!(loss, data, opt)
end

# Let's check the final values of the parameters:

params(model)

# ## Visualizing the results

#-

# What does this neural network represent? It is simply a more complicated function with two inputs and three outputs, i.e. a function $f: \mathbb{R}^2 \to \mathbb{R}^3$.
# Before, with a single layer, each component of the function $f$ basically corresponded to a hyperplane; now it will instead be a **more complicated nonlinear function** of the input data!

#-

# #### Exercise 3
#
# Visualize each component of the output separately as a heatmap and/or contours superimposed on the data. Interpret the results.

#-

# #### Solution

coords = 0:0.01:0.8

heatmap(coords, coords, (x,y)->model([x,y]).data[1])
#contour!(coords, coords, (x,y)->model([x,y]).data[1], levels=[0.5, 0.501], lw=3)

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas")
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes")

xlims!(0.4, 0.8)
ylims!(0.1, 0.5)

# We see that the first component, which is supposed to separate apples from non-apples, has been able to learn a set that has a much more complicated shape than simply a hyperplane: the hyperplane has been bent round. Sometimes it's able to encapsulate all the apple data, while sometimes it isn't, depending on how well it learnt.

coords = 0:0.01:1

#contour(coords, coords, (x,y)->model([x,y]).data[2], levels=[0.5, 0.501], lw=3)
heatmap(coords, coords, (x,y)->model([x,y]).data[2])

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas")
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes")

xlims!(0.4, 0.8)
ylims!(0.1, 0.5)

# The second component is very successful at separating out the bananas. Since there are some apples mixed in there, it can't be expected to do too much better.

coords = 0:0.01:1

heatmap(coords, coords, (x,y)->model([x,y]).data[3])
contour!(coords, coords, (x,y)->model([x,y]).data[3], levels=[0.5, 0.501], lw=3)

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas")
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes")

xlims!(0.4, 0.8)
ylims!(0.1, 0.5)

# The third component separates grapes from the rest pretty successfully.

#-

# ## What we have learned

#-

# Adding an intermediate layer allows the network to start to deform the separating surfaces that it is learning into more complicated, nonlinear (curved) shapes. This allows it to separate data that were previously unable to be separated!
#
# However, using only two features means that data from different classes overlaps. To distinguish it we would need to use more features.

#-

# ### Exercise 4
#
# Use three features (red, green and blue) and build a network with one hidden layer. Does this help to distinguish the data better?

max_stars_repo_path: Courses/Deep learning with Flux/1700.Introduction-to-deep-learning-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


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

using Markdown
using InteractiveUtils

# ╔═╡ 6a2c9516-dbda-4005-8d85-c19b5c0f54ee
begin
  using Pkg
  Pkg.activate(Base.Filesystem.homedir() * "/.config/julia/projects/oft")
  using Plots
  using LinearAlgebra
end

# ╔═╡ 4cdf2cb0-568c-47f6-bb13-a4a693e042f6
md"""
## Affine Transformations
An _affine transformation_ is a _linear transformation_ plus _a displacment_.
```math
\begin{array}
  &\mathbb{R}^2 &\to &\mathbb{R}^2& \\
  \begin{bmatrix}
    t_{1}' \\
    t_{2}'
  \end{bmatrix}
  &\mapsto
  &A
  \begin{bmatrix}
    t_{1} \\
    t_{2}
\end{bmatrix}
  + b&
  =
  \begin{bmatrix}
    a_{1,1} & a_{1,2} \\
    a_{2,1} & a_{2,2} \\
  \end{bmatrix}
  \begin{bmatrix}
    t_{1} \\
    t_{2}
  \end{bmatrix}
  +
  \begin{bmatrix}
    b_{1} \\
    b_{2}
  \end{bmatrix}
\end{array}
```
"""

# ╔═╡ a7ce1316-bac0-47cc-9829-c8d3062d3169
# begin
#   flag_x = [-1, -1, 7, 4, 7, -1]
#   flag_y = [-5, 5, 5, 3, 1, 1]
#   flag = plot(flag_x,
# 		      flag_y,
# 		      lw=3,
# 		      framestyle=:zerolines,
# 	          aspect_ratio=:equal,
# 		      linecolor=:red,
# 		      legend=false,
# 	)
# end

# ╔═╡ 019427cd-dee6-4d81-9e0e-ba86657d2d12
begin
  #flag_x = [-2, -2, 5, 3, 5, -2]
  #flag_y = [-5, 5, 5, 3, 1, 1]
  flag_xy = transpose([-2 -2 5 3 5 -2
                       -5  5 5 3 1  1])
  flag = plot(flag_xy[1:end, 1],
              flag_xy[1:end, 2],
              lw=3,
              framestyle=:zerolines,
              aspect_ratio=:equal,
              linecolor=:red,
              legend=false,
              xlims=[-10,10],
              ylims=[-10,10],
  )
end

# ╔═╡ c293aa81-9192-4cc2-8221-55077399f475
#displacement = [2 2]
displacement = [[4] [3]]

# ╔═╡ 6819c041-ad77-4ede-bf8e-32d67b589fa8
flag

# ╔═╡ ecf18f1d-7f6e-486e-b434-529787f5ef3f
begin
  ## No use.
  #flag

  plot!(flag_xy[1:end, 1],
       flag_xy[1:end, 2],
              lw=3,
              framestyle=:zerolines,
              aspect_ratio=:equal,
              linecolor=:gray,
              legend=false,
              xlims=[-10,10],
              ylims=[-10,10],
  )
  flag_displaced = plot!((flag_xy .+ displacement)[1:end, 1],
                        (flag_xy .+ displacement)[1:end, 2],
                        lw=3,
                        framestyle=:zerolines,
                        aspect_ratio=:equal,
                        linecolor=:red,
                        legend=false,
                        xlims=[-10, 10],
                        ylims=[-10, 10],
  )
end

# ╔═╡ f4611f9f-7f42-4164-8135-0fdcbf8f7ca0
begin
  #chapeau_xy = [-1 -1 7 4 7 -1
  #              -5  5 5 3 1  1]
  function draw_flag(A, b, linecolor)
    #plot()
    xmax = 15
    ymax = xmax
    chapeau_xy = [-1 -1 7 4 7 -1
                  -5  5 5 3 1  1]
    transformed = A*chapeau_xy .+ b
    plot!(transformed[1, 1:end],
          transformed[2, 1:end],
          lw=3,
          framestyle=:zerolines,
          aspect_ratio=:equal,
          linecolor=linecolor,
          legend=false,
          xlims=[-xmax, xmax],
          ylims=[-ymax, ymax],
    )
  end
end

# ╔═╡ 9c16f909-cfa8-4230-ae2c-65a599f7265d
I(3)

# ╔═╡ 58e1a38d-e958-43ff-b5f1-1a80dbe9b840
I(2) * [-1 -1 7 4 7 -1
 -5  5 5 3 1  1]

# ╔═╡ cf03b357-bd2a-4f41-a7b0-aa31628b1df6
[-1 -1 7 4 7 -1
 -5  5 5 3 1  1] .+ [4 5]'

# ╔═╡ 6554aaa3-50df-4847-b568-75a2ff5852dc
typeof([4; 5])

# ╔═╡ 5281457f-b58c-464b-97a6-efc50de31443
typeof([4, 5])

# ╔═╡ ca8edc5c-f713-4584-99f9-3fd843d31b95
typeof([4 5])

# ╔═╡ 7145ab5c-3aff-4767-9fa6-46de9f46e29d
[4 5]

# ╔═╡ adbf610d-7cb5-49c6-92b3-ca1ee519e51a
[4 5]'

# ╔═╡ 6aa15192-3ec7-417f-9847-97248d7b3ca5
draw_flag(I(2), [4 3]', :red)

# ╔═╡ a8f0d2e2-59fa-46af-bafe-5b67b799a1f6
I(2) * [-1 -1 7 4 7 -1
        -5  5 5 3 1  1] .+ [3 4]'

# ╔═╡ Cell order:
# ╠═6a2c9516-dbda-4005-8d85-c19b5c0f54ee
# ╟─4cdf2cb0-568c-47f6-bb13-a4a693e042f6
# ╠═a7ce1316-bac0-47cc-9829-c8d3062d3169
# ╠═019427cd-dee6-4d81-9e0e-ba86657d2d12
# ╠═c293aa81-9192-4cc2-8221-55077399f475
# ╠═6819c041-ad77-4ede-bf8e-32d67b589fa8
# ╠═ecf18f1d-7f6e-486e-b434-529787f5ef3f
# ╠═f4611f9f-7f42-4164-8135-0fdcbf8f7ca0
# ╠═9c16f909-cfa8-4230-ae2c-65a599f7265d
# ╠═58e1a38d-e958-43ff-b5f1-1a80dbe9b840
# ╠═cf03b357-bd2a-4f41-a7b0-aa31628b1df6
# ╠═6554aaa3-50df-4847-b568-75a2ff5852dc
# ╠═5281457f-b58c-464b-97a6-efc50de31443
# ╠═ca8edc5c-f713-4584-99f9-3fd843d31b95
# ╠═7145ab5c-3aff-4767-9fa6-46de9f46e29d
# ╠═adbf610d-7cb5-49c6-92b3-ca1ee519e51a
# ╠═6aa15192-3ec7-417f-9847-97248d7b3ca5
# ╠═a8f0d2e2-59fa-46af-bafe-5b67b799a1f6


text:
module GraphTools

#__precompile__(true)

export @graph
export +
export draw
export GraphStore

include("graphstore.jl")

end