text:
n = 1000
targets = rand(0:1, n)
scores = rand(n)
thres = 0.7
predicts = scores .>= thres

p  = sum(targets .== 1)
n  = sum(targets .== 0)
tp = sum(targets .* predicts)
tn = sum((1 .- targets) .* (1 .- predicts))
fp = sum((1 .- targets) .* predicts)
fn = sum(targets .* (1 .- predicts))

cm = ConfusionMatrix(p, n, tp, tn, fp, fn)
cm2 = ConfusionMatrix(2*p, 2*n, 2*tp, 2*tn, 2*fp, 2*fn)

set_encoding(OneZero())

encs = [
    OneZero(),
    OneMinusOne(),
    OneTwo(),
    OneVsOne(3,4),
    OneVsOne(:three,:four),
    OneVsOne("three","four"),
    OneVsRest(1, [2,3,4]),
    OneVsRest(:one, [:two, :three]),
    OneVsRest("one", ["two", "three"]),
    RestVsOne([1,2,3], 4),
    RestVsOne([:one, :two], :three),
    RestVsOne(["one", "two"], "three")
]

@testset "ConfusionMatrix +" begin
    @test cm + cm == cm2
end

@testset "ConfusionMatrix constructors for $enc encoding" for enc in encs
    global targets, predicts

    targets = recode.(current_encoding(), enc, targets)
    predicts = recode.(current_encoding(), enc, predicts)
    set_encoding(enc)

    @test ConfusionMatrix(targets, predicts) == cm
    @test ConfusionMatrix(enc, targets, predicts) == cm
    @test ConfusionMatrix(targets, scores, thres) == cm
    @test ConfusionMatrix(enc, targets, scores, thres) == cm
    @test ConfusionMatrix(targets, scores, [thres]) == [cm]
    @test ConfusionMatrix(enc, targets, scores, [thres]) == [cm]
end


text:
module SolidModeling

using LinearAlgebra

const PlaneEpsilon = 1e-8

include("structs.jl")

function fromPolygons(polygons)
    csg = Solid(Array{Vertex,1}(), Array{Int64,1}())

    p::Int64 = 0
    for i = 0:(length(polygons) - 1)
        poly = polygons[i + 1]
        for j = 2:(length(poly.vertices) - 1)
            push!(csg.vertices, poly.vertices[1])
            push!(csg.indices, p)
            push!(csg.vertices, poly.vertices[j])
            push!(csg.indices, p + 1)
            push!(csg.vertices, poly.vertices[j + 1])
            push!(csg.indices, p + 2)

            p = p + 3
        end
    end
    return csg
end

function toPolygons(model)
    list = Array{Polygon,1}();
    for i = 0:3:(length(model.indices) - 1)

        triangle = Array{Vertex,1}()
        for j = 0:2
            v = model.vertices[model.indices[i + j + 1] + 1]
            push!(triangle, copy(v));
        end
        push!(list, Polygon(triangle, fromPoints(triangle[1].pos, triangle[2].pos, triangle[3].pos)))
    end
    return list;
end

"""
    bunion(first::Solid, second::Solid)

Return new solid after union of `A` and `B`.

    +-------+            +-------+
    |       |            |       |
    |   A   |            |       |
    |    +--+----+   =   |       +----+
    +----+--+    |       +----+       |
         |   B   |            |       |
         |       |            |       |
         +-------+            +-------+

# Examples
```julia
c1 = cube(0.5, 0.5, 0.5, 1.0, 1.0, 1.0)
c2 = cube(0.5, 0.5, 0.5, 2.0, 2.0, 2.0)
r = bunion(c1, c2) # union of solids c1 and c2

v = volume(r) # volume of union c1 c2
```
"""
function bunion(first::Solid, second::Solid)
    a = Node(nothing, nothing, nothing, Array{Polygon,1}())
    b = Node(nothing, nothing, nothing, Array{Polygon,1}())

    build(a, toPolygons(first));
    build(b, toPolygons(second));

    clipTo(a, b)
    clipTo(b, a)
    invert(b)
    clipTo(b, a)
    invert(b)

    build(a, allPolygons(b))
    return fromPolygons(allPolygons(a))
end

"""
    bsubtract(first::Solid, second::Solid)

Return new solid after performing `B`-`A`.

    +-------+            +-------+
    |       |            |       |
    |   A   |            |       |
    |    +--+----+   =   |    +--+
    +----+--+    |       +----+
         |   B   |
         |       |
         +-------+

# Examples
```julia
c1 = cube(0.5, 0.5, 0.5, 1.0, 1.0, 1.0)
c2 = cube(0.5, 0.5, 0.5, 2.0, 2.0, 2.0)
r = bsubtract(c1, c2) # subtraction of c2 from c1

v = volume(r) # volume of c1-c2
```
"""
function bsubtract(first::Solid, second::Solid)
    a = Node(nothing, nothing, nothing, Array{Polygon,1}())
    b = Node(nothing, nothing, nothing, Array{Polygon,1}())

    build(a, toPolygons(first));
    build(b, toPolygons(second));

    invert(a)
    clipTo(a, b)
    clipTo(b, a)
    invert(b)
    clipTo(b, a)
    invert(b)
    build(a, allPolygons(b))
    invert(a)
    return fromPolygons(allPolygons(a))
end

"""
bintersect(first::Solid, second::Solid)

Return new solid after computing intersection of `A` and `B`.

    +-------+
    |       |
    |   A   |
    |    +--+----+   =   +--+
    +----+--+    |       +--+
         |   B   |
         |       |
         +-------+

# Examples
```julia
c1 = cube(0.5, 0.5, 0.5, 1.0, 1.0, 1.0)
c2 = cube(0.5, 0.5, 0.5, 2.0, 2.0, 2.0)
r = bintersect(c1, c2) # intersection of c1 and c2

v = volume(r)
```
"""
function bintersect(first::Solid, second::Solid)
    a = Node(nothing, nothing, nothing, Array{Polygon,1}())
    b = Node(nothing, nothing, nothing, Array{Polygon,1}())

    build(a, toPolygons(first));
    build(b, toPolygons(second));

    invert(a)
    clipTo(b, a)
    invert(b)
    clipTo(a, b)
    clipTo(b, a)
    build(a, allPolygons(b))
    invert(a)
    return fromPolygons(allPolygons(a))
end

function interpolate(vertex, other, t)
    return Vertex([vertex.pos[1] + (other.pos[1]-vertex.pos[1])*t, vertex.pos[2] + (other.pos[2]-vertex.pos[2])*t, vertex.pos[3] + (other.pos[3]-vertex.pos[3])*t])
end

function fromPoints(a, b, c)
    cross(b - a, c - a)
    n = normalize(cross(b - a, c - a))
    return Plane(n, dot(n, a))
end

function splitPolygon(plane, polygon, coplanarFront, coplanarBack, front, back)
    COPLANAR = 0
    FRONT = 1
    BACK = 2
    SPANNING = 3

    # Classify each point as well as the entire polygon into one of the above four classes.
    polygonType = 0
    types = Array{Int64,1}()

    for vertex in polygon.vertices
        t = dot(plane.normal, vertex.pos) - plane.w
        type::Int64 = (t < -PlaneEpsilon) ? BACK : ((t > PlaneEpsilon) ? FRONT : COPLANAR)
        polygonType |= type
        push!(types, type)
    end

    if polygonType === COPLANAR
        push!(dot(plane.normal, polygon.plane.normal) > 0 ? coplanarFront : coplanarBack, polygon)
    elseif polygonType === FRONT
        push!(front, polygon)
    elseif polygonType === BACK
        push!(back, polygon)
    elseif polygonType === SPANNING
        f = Array{Vertex,1}()
        b = Array{Vertex,1}()

        for i = 0:(length(polygon.vertices) - 1)
            j = (i + 1) % length(polygon.vertices)
            ti = types[i + 1]
            tj = types[j + 1]
            vi = polygon.vertices[i + 1]
            vj = polygon.vertices[j + 1]

            if ti != BACK push!(f, copy(vi)) end
            if ti != FRONT push!(b, copy(vi)) end
            if (ti | tj) == SPANNING
                t = (plane.w - dot(plane.normal, vi.pos)) / dot(plane.normal, vj.pos - vi.pos)
                v = interpolate(vi, vj, t)
                push!(f, v)
                push!(b, v)
            end
        end

        if (length(f) >= 3) push!(front, Polygon(f, fromPoints(f[1].pos, f[2].pos, f[3].pos))) end
        if (length(b) >= 3) push!(back, Polygon(b, fromPoints(b[1].pos, b[2].pos, b[3].pos))) end
    end
end

function invert(node::Node)
    for poly in node.polygons
        reverse!(poly.vertices)
        poly.plane.normal = -1 * poly.plane.normal
        poly.plane.w = -1 * poly.plane.w
    end

    if node.plane !== nothing
        node.plane.normal = -1 * node.plane.normal
        node.plane.w = -1 * node.plane.w
    end

    if node.front !== nothing invert(node.front) end
    if node.back !== nothing invert(node.back) end

    # swap front and back in place
    node.front, node.back = node.back, node.front
end

function clipPolygons(node::Node, polygons::Array{Polygon,1})
    polysize = length(polygons)

    if node.plane === nothing return copy(polygons) end

    front = Array{Polygon,1}()
    back = Array{Polygon,1}()

    for polygon in polygons
        splitPolygon(node.plane, polygon, front, back, front, back);
    end

    if node.front !== nothing
        front = clipPolygons(node.front, front)
    end

    if node.back !== nothing
        back = clipPolygons(node.back, back)
    else
        back = Array{Polygon,1}()
    end

    return vcat(front, back)
end

function clipTo(node::Node, bsp)
    node.polygons = clipPolygons(bsp, node.polygons)

    if node.front !== nothing clipTo(node.front, bsp) end
    if node.back !== nothing clipTo(node.back, bsp) end
end

function allPolygons(node::Node)
    polygons = deepcopy(node.polygons)

    if node.front !== nothing polygons = vcat(polygons, allPolygons(node.front)) end
    if node.back !== nothing polygons = vcat(polygons, allPolygons(node.back)) end

    return polygons
end

function build(node::Node, polygons)
    if length(polygons) === 0 return end

    if node.plane === nothing node.plane = deepcopy(polygons[1].plane) end

    front = Array{Polygon,1}()
    back = Array{Polygon,1}()

    for polygon in polygons
        splitPolygon(node.plane, polygon, node.polygons, node.polygons, front, back)
    end

    if length(front) > 0
        if node.front === nothing node.front = Node(nothing, nothing, nothing, Array{Polygon,1}()) end
        build(node.front, front)
    end

    if length(back) > 0
        if node.back === nothing node.back = Node(nothing, nothing, nothing, Array{Polygon,1}()) end
        build(node.back, back)
    end
end

function cube(center::Vector{Float64}, lx::Float64, ly::Float64, lz::Float64)::Solid
    return cube(center[1] - lx / 2, center[2] - ly / 2, center[3] - lz / 2, center[1] + lx / 2, center[2] + ly / 2, center[3] + lz / 2)
end

"""
    cube(xMin::Float64, yMin::Float64, zMin::Float64, xMax::Float64, yMax::Float64, zMax::Float64)::Solid
    cube(center::Vector{Float64}, lx::Float64, ly::Float64, lz::Float64)::Solid

Creates a cube by specifying its bounding box `xMin`, `yMin`, `zMin` and `xMax`, `yMax`, `zMax`.

Or create a cube using a center point located in `center` and dimensions `lx`, `ly`, and `lz`.

# Examples
```julia-repl
julia> cube(0, 0, 0, 1, 1, 1)
Main.SolidModeling.Solid(...)

julia> cube([0.5, 0.5, 0.5], 1, 1, 1)
Main.SolidModeling.Solid(...)
```
"""
function cube(xMin::Float64, yMin::Float64, zMin::Float64, xMax::Float64, yMax::Float64, zMax::Float64)::Solid
    v1 = Vertex([xMin, yMin, zMax])
    v2 = Vertex([xMin, yMax, zMax])
    v3 = Vertex([xMax, yMax, zMax])
    v4 = Vertex([xMax, yMin, zMax])
    v5 = Vertex([xMin, yMin, zMin])
    v6 = Vertex([xMin, yMax, zMin])
    v7 = Vertex([xMax, yMax, zMin])
    v8 = Vertex([xMax, yMin, zMin])

    f1p = Polygon([v1, v4, v3, v2], fromPoints(v1.pos, v4.pos, v3.pos))
    f2p = Polygon([v7, v8, v5, v6], fromPoints(v7.pos, v8.pos, v5.pos))
    f3p = Polygon([v1, v2, v6, v5], fromPoints(v1.pos, v2.pos, v6.pos))
    f4p = Polygon([v2, v3, v7, v6], fromPoints(v2.pos, v3.pos, v7.pos))
    f5p = Polygon([v3, v4, v8, v7], fromPoints(v3.pos, v4.pos, v8.pos))
    f6p = Polygon([v4, v1, v5, v8], fromPoints(v4.pos, v1.pos, v5.pos))

    polys = [f1p, f2p, f3p, f4p, f5p, f6p];

    return fromPolygons(polys);
end

function signedVolumeOfTriangle(p1::Array{Float64}, p2::Array{Float64}, p3::Array{Float64})
   	v321 = p3[1] * p2[2] * p1[3]
   	v231 = p2[1] * p3[2] * p1[3]
   	v312 = p3[1] * p1[2] * p2[3]
   	v132 = p1[1] * p3[2] * p2[3]
    v213 = p2[1] * p1[2] * p3[3]
    v123 = p1[1] * p2[2] * p3[3]

   	return (1.0 / 6.0) * (-v321 + v231 + v312 - v132 - v213 + v123)
end

"""
    volume(c::Solid)::Float64

Return the volume of **Solid** `c`.

# Examples
```julia-repl
julia> volume(cube(0.0, 0.0, 0.0, 1.0, 1.0, 1.0))
0.9999999999999999
```
"""
function volume(c::Solid)::Float64
    volume = 0.0;

    for i = 0:3:(length(c.indices) - 1)
        dv = signedVolumeOfTriangle(c.vertices[i + 1].pos, c.vertices[i + 2].pos, c.vertices[i + 3].pos)
        volume += dv
    end

    return volume
end

export cube;
export bunion;
export bsubtract;
export bintersect;
export volume;

end # module

ext: jl


text:
using DifferentialForms

using LinearAlgebra
using StaticArrays
using Test

@testset "Multivector masks and indices D=$D" for D in 0:Dmax
    Mbits = 2^D
    Mmax = UInt64(2)^Mbits - 1
    for iter in 1:100
        M = rand(1:Mmax)        # Don't test M == 0
        for iter2 in 1:100
            lin = rand(1:Multivectors.numelts(Val(D), Val(M)))
            bits = Multivectors.lin2bit(Val(D), Val(M), lin)
            @test Multivectors.bit2lin(Val(D), Val(M), bits) == lin
            lst = Multivectors.bit2lst(Val(D), bits)
            @test Multivectors.lst2bit(Val(D), lst) == bits
            @test Multivectors.lin2lst(Val(D), Val(M), lin) == lst
            @test Multivectors.lst2lin(Val(D), Val(M), lst) == lin
        end
    end
end

@testset "Create multivectors D=$D" for D in 0:Dmax
    Mbits = 2^D
    Mmax = UInt64(2)^Mbits - 1

    T = BigRat

    for iter in 1:(100 ÷ (D + 1))
        M = rand(0:Mmax)
        N = length(Multivector{D,M})
        x = Multivector{D,M,T}(SVector{N,T}(rand(T) for n in 1:N))
        X = typeof(x)
        x′ = X(SVector{N,T}(rand(T) for n in 1:N))

        @test Multivector{D,M,Float64}(x) == Multivector{D,M,Float64}(x.elts)

        @test Multivector{D,M,Float64}(collect(x)) == Multivector{D,M,Float64}(x)
    end
end

@testset "Vector space operations on multivectors D=$D" for D in 0:Dmax
    Mbits = 2^D
    Mmax = UInt64(2)^Mbits - 1

    T = Rational{Int64}

    for iter in 1:(100 ÷ (D + 1))
        Mn = zero(Mmax)
        Mx = rand(0:Mmax)
        My = rand(0:Mmax)
        Mz = rand(0:Mmax)
        n = zero(Multivector{D,Mn,T})
        x = rand(Multivector{D,Mx,T})
        y = rand(Multivector{D,My,T})
        z = rand(Multivector{D,Mz,T})
        a = rand(T)
        b = rand(T)

        # Vector space

        @test +x == x
        @test (x + y) + z == x + (y + z)

        @test n + x == x
        @test x + n == x

        @test x + y == y + x

        @test x + (-x) == n
        @test (-x) + x == n
        @test x - y == x + (-y)

        @test (a * b) * x == a * (b * x)
        @test x * (a * b) == (x * a) * b

        @test one(T) * a == a
        @test one(T) * x == x
        @test x * one(T) == x

        @test zero(T) * a == zero(T)
        @test zero(T) * x == n
        @test x * zero(T) == n

        if a != zero(T)
            @test a * inv(a) == one(T)

            @test inv(a) * (a * x) == x
            @test (x * a) * inv(a) == x

            @test inv(a) * x == a \ x
            @test x * inv(a) == x / a
        end

        @test (a + b) * x == a * x + b * x
        @test x * (a + b) == x * a + x * b

        @test a * (x + y) == a * x + a * y
        @test (x + y) * a == x * a + y * a

        # Nonlinear transformations
        @test map((a -> 2a), x) == 2 * x
        @test map((a, b) -> 2a + b, x, -x) == x
        @test map((a, b, c) -> 2a + b - c, x, x, 2 * x) == x
    end
end

@testset "Multivector algebra D=$D" for D in 0:Dmax
    Mbits = 2^D
    Mmax = UInt64(2)^Mbits - 1

    T = Rational{Int64}

    for iter in 1:(100 ÷ (D + 1))
        Me = UInt64(0b1)
        Mx = rand(0:Mmax)
        My = rand(0:Mmax)
        Mz = rand(0:Mmax)
        e = one(Multivector{D,Me,T})
        x = rand(Multivector{D,Mx,T})
        y = rand(Multivector{D,My,T})
        z = rand(Multivector{D,Mz,T})
        a = rand(T)
        b = rand(T)

        # Multiplicative structure

        # units
        if D == 2
            e = unit(Multivector{D,Me,T})
            e1 = unit(Multivector{D,UInt64(0b0010),T}, 1)
            e2 = unit(Multivector{D,UInt64(0b00100),T}, 2)
            e12 = unit(Multivector{D,UInt64(0b1000),T}, 1, 2)

            @test ~e == e
            @test ~e1 == e1
            @test ~e2 == e2
            @test ~e12 == -e12

            @test ⋆e == e12
            @test ⋆e1 == e2
            @test ⋆e2 == -e1
            @test ⋆e12 == e

            @test e ∧ e == e
            @test e ∧ e1 == e1
            @test e ∧ e2 == e2
            @test e ∧ e12 == e12
            @test e1 ∧ e == e1
            @test e2 ∧ e == e2
            @test e12 ∧ e == e12
            @test e1 ∧ e1 == 0 * e12
            @test e1 ∧ e2 == e12
            @test e2 ∧ e1 == -e12
            @test e2 ∧ e2 == 0 * e12

            @test e ∨ e12 == e
            @test e1 ∨ e1 == 0 * e
            @test e2 ∨ e2 == 0 * e
            @test e1 ∨ e2 == e
            @test e1 ∨ e12 == e1
            @test e2 ∨ e12 == e2
            @test e12 ∨ e12 == e12

            @test e ⋅ e == e
            @test e1 ⋅ e1 == e
            @test e2 ⋅ e2 == e
            @test e1 ⋅ e2 == 0 * e
            @test e12 ⋅ e12 == e

            @test e × e == e12
            @test e × e1 == e2
            @test e × e2 == -e1
            @test e × e12 == e
            @test e1 × e1 == 0 * e
            @test e2 × e2 == 0 * e
            @test e1 × e2 == e

            @test e * e == e
            @test e * e1 == e1
            @test e * e2 == e2
            @test e * e12 == e12
            @test e1 * e == e1
            @test e2 * e == e2
            @test e12 * e == e12
            @test e1 * e1 == e
            @test e1 * e2 == e12
            @test e2 * e1 == -e12
            @test e2 * e2 == e
        end

        #TODO     # various duals
        #TODO     cycle_basis(x)::typeof(x)
        #TODO     x′ = x
        #TODO     for d in 1:D
        #TODO         x′ = cycle_basis(x′)
        #TODO     end
        #TODO     @test x′ == x
        #TODO     reverse_basis(x)::typeof(x)
        #TODO     @test reverse_basis(reverse_basis(x)) == x
        #TODO     @test reverse_basis(cycle_basis(reverse_basis(cycle_basis(x)))) == x

        @test ~~x == x
        @test a * ~x == ~(a * x)

        @test conj(conj(x)) == x
        @test a * conj(x) == conj(a * x)

        # @test ⋆⋆x == bitsign(R1 * (D - R1)) * x
        @test inv(⋆)(⋆x) == x
        @test ⋆inv(⋆)(x) == x
        @test ⋆⋆ ⋆ ⋆x == x
        @test a * ⋆x == ⋆(a * x)

        # exterior product: x ∧ y
        @test ∧(x) == x

        (x ∧ y)::Multivector{D}

        (x ∧ y ∧ z)::Multivector{D}
        @test (x ∧ y) ∧ z == x ∧ (y ∧ z)
        @test x ∧ y ∧ z == (x ∧ y) ∧ z
        @test ∧(x, y, z) == (x ∧ y) ∧ z
        @test ∧(SVector{0,typeof(x)}()) == e
        @test ∧(SVector(x)) == x
        @test ∧(SVector(x, y)) == x ∧ y
        @test ∧(SVector(x, y, z)) == (x ∧ y) ∧ z
        @test ∧(SVector{0,typeof(x)}()) ∧ x == x
        @test x ∧ ∧(SVector{0,typeof(x)}()) == x

        @test e ∧ x == x
        @test x ∧ e == x

        @test x ∧ zero(y) == zero(x ∧ y)
        @test zero(y) ∧ x == zero(y ∧ x)

        @test a * (x ∧ y) == x ∧ (a * y)
        @test x ∧ (y + z) == x ∧ y + x ∧ z

        @test x ∧ y == ~(~y ∧ ~x)

        # regressive product: ⋆(x ∨ y) = ⋆x ∧ ⋆y
        @test ∨(x) == x
        (x ∨ y)::Multivector{D}
        @test ⋆(x ∨ y) == ⋆x ∧ ⋆y

        (x ∨ y ∨ z)::Multivector{D}
        @test ⋆(x ∨ y ∨ z) == ⋆x ∧ ⋆y ∧ ⋆z
        @test x ∨ y ∨ z == (x ∨ y) ∨ z
        @test ∨(x, y, z) == (x ∨ y) ∨ z
        @test ∨(SVector{0,typeof(x)}()) == ⋆e
        @test ∨(SVector(x)) == x
        @test ∨(SVector(x, y)) == x ∨ y
        @test ∨(SVector(x, y, z)) == (x ∨ y) ∨ z
        @test ∨(SVector{0,typeof(x)}()) ∨ x == x
        @test x ∨ ∨(SVector{0,typeof(x)}()) == x

        @test ⋆e ∨ x == x
        @test x ∨ ⋆e == x

        @test x ∨ zero(y) == zero(x ∨ y)
        @test zero(y) ∨ x == zero(y ∨ x)

        @test a * (x ∨ y) == x ∨ (a * y)
        @test x ∨ (y + z) == (x ∨ y) + (x ∨ z)

        Q = inv(⋆) ∘ (~) ∘ (⋆)
        Q1 = (⋆) ∘ (~) ∘ inv(⋆)
        @test x ∨ y == Q1(Q(y) ∨ Q(x))

        # dot product: x ⋅ y = x ∨ ⋆y   (right contraction)
        (x ⋅ y)::Multivector{D}
        @test x ⋅ y == x ∨ ⋆y

        norm2(x)::T
        @test norm2(x) >= 0
        @test norm2(a * x) == norm2(a) * norm2(x)
        # @test norm2(y + y2) <= norm2(y) + norm2(y2)
        norm(x)::float(T)
        @test norm(x) ≈ sqrt(norm2(x))
        @test norm(a * x) ≈ norm(a) * norm(x)
        @test norm(x + y) <= norm(x) + norm(y) || norm(x + y) ≈ norm(x) + norm(y)

        # cross product: x × y = ⋆(x ∧ y)
        (x × y)::Multivector{D}
        @test x × y == ⋆(x ∧ y)

        # geometric product: x * y
        @test *(x) == x

        (x * y)::Multivector{D}

        (x * y * z)::Multivector{D}
        @test (x * y) * z == x * (y * z)
        @test x * y * z == (x * y) * z
        @test *(x, y, z) == (x * y) * z
        @test *(SVector{0,typeof(x)}()) == e
        @test *(SVector(x)) == x
        @test *(SVector(x, y)) == x * y
        @test *(SVector(x, y, z)) == (x * y) * z
        @test *(SVector{0,typeof(x)}()) * x == x
        @test x * *(SVector{0,typeof(x)}()) == x

        @test e * x == x
        @test x * e == x

        @test x * zero(y) == zero(x * y)
        @test zero(y) * x == zero(y * x)

        @test a * (x * y) == x * (a * y)
        @test x * (y + z) == x * y + x * z
    end
end

#TODO const Dmax4 = min(4, Dmax)
#TODO @testset "Tensor sums of multivectors D1=$D1 D2=$D2 R=$R" for D1 in 0:Dmax4,
#TODO D2 in 0:Dmax4,
#TODO R in 1:min(D1, D2)
#TODO
#TODO     D = D1 + D2
#TODO
#TODO     T = Rational{Int64}
#TODO     x = rand(Form{D1,R,T})
#TODO     x2 = rand(Form{D1,R,T})
#TODO     y = rand(Form{D2,R,T})
#TODO     y2 = rand(Form{D2,R,T})
#TODO     a = rand(T)
#TODO
#TODO     # units
#TODO     if D1 == 1 && D2 == 1 && R == 1
#TODO         u(d, inds...) = unit(Form{d,length(inds),T}, inds...)
#TODO
#TODO         @test u(1, 1) ⊕ u(1, 1) == u(2, 1) + u(2, 2)
#TODO
#TODO         @test u(1, 1) ⊕ u(2, 1) == u(3, 1) + u(3, 2)
#TODO         @test u(1, 1) ⊕ u(2, 2) == u(3, 1) + u(3, 3)
#TODO
#TODO         @test u(3, 1) ⊕ u(2, 1) == u(5, 1) + u(5, 4)
#TODO         @test u(3, 1) ⊕ u(2, 2) == u(5, 1) + u(5, 5)
#TODO         @test u(3, 2) ⊕ u(2, 1) == u(5, 2) + u(5, 4)
#TODO         @test u(3, 2) ⊕ u(2, 2) == u(5, 2) + u(5, 5)
#TODO         @test u(3, 3) ⊕ u(2, 1) == u(5, 3) + u(5, 4)
#TODO         @test u(3, 3) ⊕ u(2, 2) == u(5, 3) + u(5, 5)
#TODO     end
#TODO
#TODO     @test ⊕(x) == x
#TODO     @test a * ⊕(x) == ⊕(a * x)
#TODO
#TODO     (x ⊕ y)::Form{D,R,T}
#TODO     @test a * (x ⊕ y) == (a * x) ⊕ (a * y)
#TODO     @test (x ⊕ y) * a == (x * a) ⊕ (y * a)
#TODO     @test x ⊕ y == (zero(x) ⊕ y) + (x ⊕ zero(y))
#TODO     @test (x + x2) ⊕ (y + y2) == (x ⊕ y) + (x2 ⊕ y2)
#TODO
#TODO     @test reverse_basis(x ⊕ y) ==
#TODO           # bitsign(R * R) *
#TODO           (reverse_basis(y) ⊕ reverse_basis(x))
#TODO
#TODO     for D3 in R:Dmax4
#TODO         z = rand(Form{D3,R,T})
#TODO         @test (x ⊕ y) ⊕ z == x ⊕ (y ⊕ z)
#TODO         @test x ⊕ y ⊕ z == x ⊕ (y ⊕ z)
#TODO         @test ⊕(x, y, z) == x ⊕ (y ⊕ z)
#TODO     end
#TODO end
#TODO
#TODO const Dmax3 = min(3, Dmax)
#TODO @testset "Tensor products of multivectors D1=$D1 R1=$R1 D2=$D2 R2=$R2" for D1 in
#TODO                                                                            0:Dmax3,
#TODO D2 in 0:Dmax3,
#TODO R1 in 0:D1,
#TODO R2 in 0:D2
#TODO
#TODO     D = D1 + D2
#TODO     R = R1 + R2
#TODO
#TODO     T = Rational{Int64}
#TODO     e = one(Form{0,0,T})
#TODO     x = rand(Form{D1,R1,T})
#TODO     x2 = rand(Form{D1,R1,T})
#TODO     y = rand(Form{D2,R2,T})
#TODO     y2 = rand(Form{D2,R2,T})
#TODO     a = rand(T)
#TODO
#TODO     # units
#TODO     if D1 == 1 && D2 == 1 && R1 == 1 && R2 == 1
#TODO         u(d, inds...) = unit(Form{d,length(inds),T}, inds...)
#TODO
#TODO         @test u(1) ⊗ u(1) == u(2)
#TODO         @test u(1, 1) ⊗ u(1) == u(2, 1)
#TODO         @test u(1) ⊗ u(1, 1) == u(2, 2)
#TODO         @test u(1, 1) ⊗ u(1, 1) == u(2, 1, 2)
#TODO
#TODO         @test u(1) ⊗ u(2) == u(3)
#TODO         @test u(1, 1) ⊗ u(2) == u(3, 1)
#TODO         @test u(1) ⊗ u(2, 1) == u(3, 2)
#TODO         @test u(1, 1) ⊗ u(2, 1) == u(3, 1, 2)
#TODO         @test u(1) ⊗ u(2, 2) == u(3, 3)
#TODO         @test u(1, 1) ⊗ u(2, 2) == u(3, 1, 3)
#TODO         @test u(1) ⊗ u(2, 1, 2) == u(3, 2, 3)
#TODO         @test u(1, 1) ⊗ u(2, 1, 2) == u(3, 1, 2, 3)
#TODO     end
#TODO
#TODO     @test ⊗(x) == x
#TODO     @test a * ⊗(x) == ⊗(a * x)
#TODO
#TODO     (x ⊗ y)::Form{D,R,T}
#TODO     @test a * (x ⊗ y) == (a * x) ⊗ y
#TODO     @test (x ⊗ y) * a == x ⊗ (y * a)
#TODO     @test e ⊗ y == y
#TODO     @test x ⊗ e == x
#TODO     @test (x + x2) ⊗ y == x ⊗ y + x2 ⊗ y
#TODO     @test x ⊗ (y + y2) == x ⊗ y + x ⊗ y2
#TODO
#TODO     @test reverse_basis(x ⊗ y) ==
#TODO           bitsign(R1 * R2) * (reverse_basis(y) ⊗ reverse_basis(x))
#TODO
#TODO     for D3 in 0:Dmax3, R3 in 0:D3
#TODO         z = rand(Form{D3,R3,T})
#TODO         @test (x ⊗ y) ⊗ z == x ⊗ (y ⊗ z)
#TODO         @test x ⊗ y ⊗ z == x ⊗ (y ⊗ z)
#TODO         @test ⊗(x, y, z) == x ⊗ (y ⊗ z)
#TODO     end
#TODO end


text:
module DependencyTree

  using Boot

  include_folder(DependencyTree, @__FILE__, is_sorted=true)

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of default printing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Default" begin

    expected = """
<!DOCTYPE html>

<meta charset="UTF-8">

      Col. 1
      Col. 2
      Col. 3
      Col. 4

      1
      false
      1.0
      1

      2
      true
      2.0
      2

      3
      false
      3.0
      3

      4
      true
      4.0
      4

      5
      false
      5.0
      5

      6
      true
      6.0
      6

"""
    result = pretty_table(String, data, backend = Val(:html))
    @test result == expected

    expected = """

      Col. 1
      Col. 2
      Col. 3
      Col. 4

      1
      false
      1.0
      1

      2
      true
      2.0
      2

      3
      false
      3.0
      3

      4
      true
      4.0
      4

      5
      false
      5.0
      5

      6
      true
      6.0
      6

"""

    result = pretty_table(String, data, backend = Val(:html), standalone = false)
    @test result == expected
end

@testset "Dictionaries" begin
    dict = Dict{Int64,String}(1 => "Jan", 2 => "Feb", 3 => "Mar", 4 => "Apr",
                              5 => "May", 6 => "Jun")

    expected = """

      Keys
      Values

      Int64
      String

      1
      Jan

      2
      Feb

      3
      Mar

      4
      Apr

      5
      May

      6
      Jun

"""

    result = pretty_table(
        String,
        dict;
        backend = Val(:html),
        sortkeys = true,
        standalone = false
    )

    @test result == expected
end

@testset "Vectors" begin

    vec = 0:1:5

    expected = """

      Col. 1

      0

      1

      2

      3

      4

      5

"""

    result = pretty_table(String, vec, backend = Val(:html), standalone = false)
    @test result == expected

    expected = """

      Row
      Col. 1

      1
      0

      2
      1

      3
      2

      4
      3

      5
      4

      6
      5

"""

    result = pretty_table(
        String,
        vec;
        alignment = :c,
        backend = Val(:html),
        show_row_number = true,
        standalone = false
    )

    @test result == expected

    expected = """

      A

      B

      C

      D

      0

      1

      2

      3

      4

      5

"""

    result = pretty_table(
        String,
        vec;
        backend = Val(:html),
        header = (["A"], ["B"], ["C"], ["D"]),
        standalone = false
    )

    @test result == expected
end

@testset "Print missing, nothing, and #undef" begin

    matrix = Matrix{Any}(undef,3,3)
    matrix[1,1:2] .= missing
    matrix[2,1:2] .= nothing
    matrix[3,1]   = missing
    matrix[3,2]   = nothing

    expected = """

      Col. 1
      Col. 2
      Col. 3

      missing
      missing
      #undef

      nothing
      nothing
      #undef

      missing
      nothing
      #undef

"""

    result = pretty_table(
        String,
        matrix;
        tf = tf_html_default,
        standalone = false
    )

    @test result == expected
end

@testset "HTML escaping" begin
    header = [
        "1",
        "2"
    ]

    matrix = [
        1 "Bold"
        2 "Italics"
        3 "
Paragraph
"
    ]

    expected = """

      <span style = "color: blue;">1</span>
      <span style = "color: blue;">2</span>

      1
      <b>Bold</b>

      2
      <em>Italics</em>

      3
      <p class="myClass">Paragraph</p>

"""

    result = pretty_table(
        String,
        matrix;
        backend = Val(:html),
        header = header,
        standalone = false
    )

    @test result == expected

    expected = """

      1
      2

      1
      Bold

      2
      Italics

      3

Paragraph

"""

    result = pretty_table(
        String,
        matrix;
        allow_html_in_cells = true,
        backend = Val(:html),
        header = header,
        standalone = false
    )

    @test result == expected
end


text:
@testset "NotImplemented" begin
    @testset "NotImplemented" begin
        ni = ChainRulesCore.NotImplemented(
            @__MODULE__, LineNumberNode(@__LINE__, @__FILE__), "error"
        )
        ni2 = ChainRulesCore.NotImplemented(
            @__MODULE__, LineNumberNode(@__LINE__, @__FILE__), "error2"
        )
        x = rand()
        thunk = @thunk(x^2)

        # conjugate
        @test conj(ni) === ni

        # addition
        for a in (true, x, NoTangent(), ZeroTangent(), thunk)
            @test ni + a === ni
            @test a + ni === ni
        end
        @test +ni === ni
        @test ni + ni2 === ni
        @test ni2 + ni === ni2

        # multiplication and dot product
        @test -ni == ni
        for a in (true, x, thunk)
            @test ni * a === ni
            @test a * ni === ni
            @test dot(ni, a) === ni
            @test dot(a, ni) === ni
        end
        for a in (NoTangent(), ZeroTangent())
            @test ni * a === a
            @test a * ni === a
            @test dot(ni, a) === a
            @test dot(a, ni) === a
        end
        @test ni * ni2 === ni
        @test ni2 * ni === ni2
        @test dot(ni, ni2) === ni
        @test dot(ni2, ni) === ni2

        # broadcasting
        @test ni .* x === ni
        @test x .* ni === ni
        @test broadcastable(ni) isa Ref{typeof(ni)}

        # unsupported operations
        E = ChainRulesCore.NotImplementedException
        for a in (true, x, NoTangent(), ZeroTangent(), thunk)
            @test_throws E ni - a
            @test_throws E a - ni
        end
        @test_throws E ni - ni2
        @test_throws E ni / x
        @test_throws E x / ni
        @test_throws E ni / ni2
        @test_throws E zero(ni)
        @test_throws E zero(typeof(ni))
        @test_throws E iterate(ni)
        @test_throws E iterate(ni, nothing)
        @test_throws E adjoint(ni)
        @test_throws E transpose(ni)
        @test_throws E convert(Float64, ni)
    end

    @testset "@not_implemented" begin
        ni = @not_implemented("myerror")
        @test ni isa ChainRulesCore.NotImplemented
        @test ni.mod isa Module
        @test ni.source isa LineNumberNode
        @test ni.info == "myerror"

        info = "some info"
        ni = @not_implemented(info)
        @test ni isa ChainRulesCore.NotImplemented
        @test ni.mod isa Module
        @test ni.source isa LineNumberNode
        @test ni.info === info
    end

    @testset "NotImplementedException" begin
        ni = @not_implemented("not implemented")
        ex = ChainRulesCore.NotImplementedException(ni)
        @test ex isa ChainRulesCore.NotImplementedException
        @test ex.mod === ni.mod
        @test ex.source === ni.source
        @test ex.info === ni.info
    end

    @testset "@scalar_rule" begin
        notimplemented1(x, y) = x + y
        @scalar_rule notimplemented1(x, y) (@not_implemented("notimplemented1"), 1)

        y, ẏ = frule((NoTangent(), 1.2, 2.3), notimplemented1, 3, 2)
        @test y == 5
        @test ẏ isa ChainRulesCore.NotImplemented

        res, pb = rrule(notimplemented1, 3, 2)
        @test res == 5
        f̄, x̄1, x̄2 = pb(3.1)
        @test f̄ == NoTangent()
        @test x̄1 isa ChainRulesCore.NotImplemented
        @test x̄2 == 3.1

        notimplemented2(x, y) = (x + y, x - y)
        @scalar_rule notimplemented2(x, y) (@not_implemented("notimplemented2"), 1) (1, -1)

        y, (ẏ1, ẏ2) = frule((NoTangent(), 1.2, 2.3), notimplemented2, 3, 2)
        @test y == (5, 1)
        @test ẏ1 isa ChainRulesCore.NotImplemented
        @test ẏ2 ≈ -1.1

        res, pb = rrule(notimplemented2, 3, 2)
        @test res == (5, 1)
        f̄, x̄1, x̄2 = pb((3.1, 4.5))
        @test f̄ == NoTangent()
        @test x̄1 isa ChainRulesCore.NotImplemented
        @test x̄2 == -1.4
    end
end


text:
# we are making a Julia model if this gets called, so lets load the Julia strategy -
include("./strategy/VFFStrategy.jl")

function extract_model_name(tree_root::XMLElement)::String

    # get model, and extract name attribute -
    model_tag = get_elements_by_tagname(tree_root,"model")[1] # only one model tag -
    return attribute(model_tag,"id")
end

function make_vff_model(path_to_sbml_file::String, path_to_output_dir::String)

    # TODO: checks - are paths legit?

    # start -
    time_start = time();

    # initialize -
    vff_component_set = Array{VLProgramComponent,1}()

    @info "Parsing $(path_to_sbml_file) ... "

    # load up ETree -
    sbml_tree = parse_file(path_to_sbml_file)
    sbml_tree_root = root(sbml_tree)

    # build the filename =
    model_name = extract_model_name(sbml_tree_root) # bad code - we need to check for this method not returning what we want
    filename = "$(model_name).vff"

    @info "Starting generation of $(filename) ... "

    # build the global header -
    global_header_program_component = build_global_header_program_component(sbml_tree_root, filename)
    push!(vff_component_set, global_header_program_component)

    # build txtl section -
    txtl_program_component = build_txtl_program_component(sbml_tree_root, filename)
    push!(vff_component_set, txtl_program_component)

    # build metabolism section -
    metabolism_program_component = build_metabolism_program_component(sbml_tree_root, filename)
    push!(vff_component_set, metabolism_program_component)

    # build global header section -
    grn_program_component = build_grn_program_component(sbml_tree_root, filename)
    push!(vff_component_set, grn_program_component)

    # Dump the vff_component_set to disk -
    write_combined_program_components_to_disk("$(path_to_output_dir)", vff_component_set)

    # stop timer -
    elapsed_time = (time() - time_start)    # in ns
    @info "Done! elapsed_time: $(time() - time_start)s"

    # free the document memory -
    free(sbml_tree)
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
"""
    evolve(c::CurveProblem, solmethod=Tsit5; mdc_callback=nothing, callback=nothing, momentum_tol=1e-3,kwargs...)
Evolves a minimally disruptive curve, with curve parameters specified by curveProblem. Uses DifferentialEquations.solve() to run the ODE.
MinimallyDisruptiveCurves.jl callbacks go in the `mdc_callback` keyword argument. You can also use any DifferentialEquations.jl callbacks compatible with DifferentialEquations.solve(). They go in the `callback` keyword argument
"""
function evolve(c::CurveProblem, solmethod=Tsit5; mdc_callback=CallbackCallable[], callback=nothing, momentum_tol=1e-3,kwargs...)

    (!(eltype(mdc_callback) == CallbackCallable)) && (mdc_callback = convert(Vector{CallbackCallable}, mdc_callback))

    function merge_sols(neg, pos, p)
        t = cat(neg.t, pos.t, dims=1)
        u = cat(neg.u, pos.u, dims=1)
        return DiffEqBase.build_solution(p, solmethod, t, u)
    end

    probs = c()

    callbacks = CallbackSet(
        build_callbacks(c, mdc_callback, momentum_tol)...,
        build_callbacks(c, callback)
    )

    runn(p) = solve(p, solmethod(); callback=callbacks, kwargs...)
    (length(probs) == 1) && (sols = runn(probs[1]))

    if length(probs) == 2
        psol = Threads.@spawn runn(probs[2])
        nsol = runn(probs[1])
        nsol.u[:] = nsol.u[end:-1:1]
        nsol.t[:] = -nsol.t[end:-1:1]
        wait(psol)
        psol = psol.result
        sols = merge_sols(nsol, psol, probs[end])
    end
    return MDCSolution(sols, c.cost)
end


text:
module MPIStateArrays

using DoubleFloats
using KernelAbstractions
using LazyArrays
using LinearAlgebra
using MPI
using StaticArrays

using ..TicToc
using ..VariableTemplates: @vars, varsindex

using Base.Broadcast: Broadcasted, BroadcastStyle, ArrayStyle

# This is so we can do things like
#   similar(Array{Float64}, Int, 3, 4)
Base.similar(::Type{A}, ::Type{FT}, dims...) where {A <: Array, FT} =
    similar(Array{FT}, dims...)

using CuArrays
Base.similar(::Type{A}, ::Type{FT}, dims...) where {A <: CuArray, FT} =
    similar(CuArray{FT}, dims...)

include("CMBuffers.jl")
using .CMBuffers

cpuify(x::AbstractArray) = convert(Array, x)
cpuify(x::Real) = x

export MPIStateArray, euclidean_distance, weightedsum

"""
    MPIStateArray{FT, DATN<:AbstractArray{FT,3}, DAI1, DAV,
                  DAT2<:AbstractArray{FT,2}} <: AbstractArray{FT, 3}
"""
mutable struct MPIStateArray{
    FT,
    V,
    DATN <: AbstractArray{FT, 3},
    DAI1,
    DAV,
    Buf <: CMBuffer,
} <: AbstractArray{FT, 3}
    mpicomm::MPI.Comm
    data::DATN
    realdata::DAV

    realelems::UnitRange{Int64}
    ghostelems::UnitRange{Int64}

    vmaprecv::DAI1
    vmapsend::DAI1

    sendreq::Array{MPI.Request, 1}
    recvreq::Array{MPI.Request, 1}

    send_buffer::Buf
    recv_buffer::Buf

    nabrtorank::Array{Int64, 1}
    nabrtovmaprecv::Array{UnitRange{Int64}, 1}
    nabrtovmapsend::Array{UnitRange{Int64}, 1}

    weights::DATN

    function MPIStateArray{FT, V}(
        mpicomm,
        DA,
        Np,
        nstate,
        numelem,
        realelems,
        ghostelems,
        vmaprecv,
        vmapsend,
        nabrtorank,
        nabrtovmaprecv,
        nabrtovmapsend,
        weights;
        mpi_knows_cuda = nothing,
    ) where {FT, V}
        data = similar(DA, FT, Np, nstate, numelem)

        if isnothing(mpi_knows_cuda)
            mpi_knows_cuda = MPI.has_cuda()
        end

        if data isa Array || mpi_knows_cuda
            kind = SingleCMBuffer
        else
            kind = DoubleCMBuffer
        end
        recv_buffer = CMBuffer{FT}(DA, kind, nstate, length(vmaprecv))
        send_buffer = CMBuffer{FT}(DA, kind, nstate, length(vmapsend))

        realdata =
            view(data, ntuple(i -> Colon(), ndims(data) - 1)..., realelems)
        DAV = typeof(realdata)

        nnabr = length(nabrtorank)
        sendreq = fill(MPI.REQUEST_NULL, nnabr)
        recvreq = fill(MPI.REQUEST_NULL, nnabr)

        # If vmap is not on the device we need to copy it up (we also do not want to
        # put it up everytime, so if it's already on the device then we do not do
        # anything).
        #
        # Better way than checking the type names?
        # XXX: Use Adapt.jl vmaprecv = adapt(DA, vmaprecv)
        if typeof(vmaprecv).name != typeof(data).name
            vmaprecv =
                copyto!(similar(DA, eltype(vmaprecv), size(vmaprecv)), vmaprecv)
        end
        if typeof(vmapsend).name != typeof(data).name
            vmapsend =
                copyto!(similar(DA, eltype(vmapsend), size(vmapsend)), vmapsend)
        end
        if typeof(weights).name != typeof(data).name
            weights =
                copyto!(similar(DA, eltype(weights), size(weights)), weights)
        end

        DAI1 = typeof(vmaprecv)
        Buf = typeof(send_buffer)
        Q = new{FT, V, typeof(data), DAI1, DAV, Buf}(
            # Make sure that each MPIStateArray has its own MPI context.  This
            # allows multiple MPIStateArrays to be communicating asynchronously
            # at the same time without having to explicitly manage tags.
            MPI.Comm_dup(mpicomm),
            data,
            realdata,
            realelems,
            ghostelems,
            vmaprecv,
            vmapsend,
            sendreq,
            recvreq,
            send_buffer,
            recv_buffer,
            nabrtorank,
            nabrtovmaprecv,
            nabrtovmapsend,
            weights,
        )
        # Make sure that we have finished all outstanding data for halo
        # exchanges before the MPIStateArray is finalize.
        finalizer(Q) do x
            MPI.Waitall!(x.recvreq)
            MPI.Waitall!(x.sendreq)
        end
        return Q
    end
end

Base.fill!(Q::MPIStateArray, x) = fill!(Q.data, x)

function Base.getproperty(Q::MPIStateArray{FT, V}, sym::Symbol) where {FT, V}
    if sym ∈ V.names
        varrange = varsindex(V, sym)
        return view(realview(Q), :, varrange, :)
    else
        return getfield(Q, sym)
    end
end

"""
   MPIStateArray{FT, V}(mpicomm, DA, Np, nstate, numelem; realelems=1:numelem,
                        ghostelems=numelem:numelem-1,
                        vmaprecv=1:0,
                        vmapsend=1:0,
                        nabrtorank=Array{Int64}(undef, 0),
                        nabrtovmaprecv=Array{UnitRange{Int64}}(undef, 0),
                        nabrtovmapsend=Array{UnitRange{Int64}}(undef, 0),
                        weights)

Construct an `MPIStateArray` over the communicator `mpicomm` with `numelem`
elements, using array type `DA` with element type `FT`. `V` is an optional type
to associate the `MPIStateArray` with a given `NamedTuple` used for `Vars` (e.g,
generated by macro `@vars`). The arrays that are held in this created
`MPIStateArray` will be of size `(Np, nstate, numelem)`.

The range `realelems` is the number of elements that this mpirank owns, whereas
the range `ghostelems` is the elements that are owned by other mpiranks.
Elements are stored as 'realelems` followed by `ghostelems`.

  * `vmaprecv` is an ordered array of elements to be received from neighboring
     mpiranks.  This is a vmap index into the MPIStateArray `A[:,*,:]`.
  * `vmapsend` is an ordered array of elements to be sent to neighboring
     mpiranks.  This is a vmap index into the MPIStateArray `A[:,*,:]`.
  * `nabrtorank` is the list of neighboring mpiranks
  * `nabrtovmaprecv` is an `Array` of `UnitRange` that give the ghost data to be
    received from neighboring mpiranks (indexes into `vmaprecv`)
  * nabrtovmapsend` is an `Array` of `UnitRange` for which elements to send to
    which neighboring mpiranks indexing into the `vmapsend`
  * `weights` is an optional array which gives weight for each degree of freedom
    to be used when computing the 2-norm of the array
"""
MPIStateArray{FT}(args...; kwargs...) where {FT} =
    MPIStateArray{FT, @vars()}(args...; kwargs...)

function MPIStateArray{FT, V}(
    mpicomm,
    DA,
    Np,
    nstate,
    numelem;
    realelems = 1:numelem,
    ghostelems = numelem:(numelem - 1),
    vmaprecv = 1:0,
    vmapsend = 1:0,
    nabrtorank = Array{Int64}(undef, 0),
    nabrtovmaprecv = Array{UnitRange{Int64}}(undef, 0),
    nabrtovmapsend = Array{UnitRange{Int64}}(undef, 0),
    weights = nothing,
    mpi_knows_cuda = nothing,
) where {FT, V}

    if weights == nothing
        weights = similar(DA, FT, ntuple(j -> 0, 3))
    end
    MPIStateArray{FT, V}(
        mpicomm,
        DA,
        Np,
        nstate,
        numelem,
        realelems,
        ghostelems,
        vmaprecv,
        vmapsend,
        nabrtorank,
        nabrtovmaprecv,
        nabrtovmapsend,
        weights,
        mpi_knows_cuda = mpi_knows_cuda,
    )
end

# FIXME: should general cases be handled?
function Base.similar(
    Q::MPIStateArray{OLDFT, V},
    ::Type{A},
    ::Type{FT},
) where {A <: AbstractArray, FT <: Number, OLDFT, V}
    MPIStateArray{FT, V}(
        Q.mpicomm,
        A,
        size(Q.data)...,
        Q.realelems,
        Q.ghostelems,
        Q.vmaprecv,
        Q.vmapsend,
        Q.nabrtorank,
        Q.nabrtovmaprecv,
        Q.nabrtovmapsend,
        Q.weights,
    )
end
function Base.similar(
    Q::MPIStateArray{FT},
    ::Type{A},
) where {A <: AbstractArray, FT <: Number}
    similar(Q, A, FT)
end
function Base.similar(Q::MPIStateArray, ::Type{FT}) where {FT <: Number}
    similar(Q, typeof(Q.data), FT)
end
function Base.similar(Q::MPIStateArray{FT}) where {FT}
    similar(Q, FT)
end

Base.size(Q::MPIStateArray, x...; kw...) = size(Q.realdata, x...; kw...)

Base.getindex(Q::MPIStateArray, x...; kw...) = getindex(Q.realdata, x...; kw...)

Base.setindex!(Q::MPIStateArray, x...; kw...) =
    setindex!(Q.realdata, x...; kw...)

Base.eltype(Q::MPIStateArray, x...; kw...) = eltype(Q.data, x...; kw...)

Base.Array(Q::MPIStateArray) = Array(Q.data)

# broadcasting stuff

# find the first MPIStateArray among `bc` arguments
# based on https://docs.julialang.org/en/v1/manual/interfaces/#Selecting-an-appropriate-output-array-1
find_mpisa(bc::Broadcasted) = find_mpisa(bc.args)
find_mpisa(args::Tuple) = find_mpisa(find_mpisa(args[1]), Base.tail(args))
find_mpisa(x) = x
find_mpisa(a::MPIStateArray, rest) = a
find_mpisa(::Any, rest) = find_mpisa(rest)

Base.BroadcastStyle(::Type{<:MPIStateArray}) = ArrayStyle{MPIStateArray}()
function Base.similar(
    bc::Broadcasted{ArrayStyle{MPIStateArray}},
    ::Type{FT},
) where {FT}
    similar(find_mpisa(bc), FT)
end

# transform all arguments of `bc` from MPIStateArrays to Arrays
function transform_broadcasted(bc::Broadcasted, ::Array)
    transform_array(bc)
end
function transform_array(bc::Broadcasted)
    Broadcasted(bc.f, transform_array.(bc.args), bc.axes)
end
transform_array(mpisa::MPIStateArray) = mpisa.realdata
transform_array(x) = x

Base.copyto!(dest::Array, src::MPIStateArray) = copyto!(dest, src.data)

function Base.copyto!(dest::MPIStateArray, src::MPIStateArray)
    copyto!(dest.realdata, src.realdata)
    dest
end

@inline function Base.copyto!(dest::MPIStateArray, bc::Broadcasted{Nothing})
    # check for the case a .= b, where b is an array
    if bc.f === identity && bc.args isa Tuple{AbstractArray}
        if bc.args isa Tuple{MPIStateArray}
            realindices = CartesianIndices((
                axes(dest.data)[1:(end - 1)]...,
                dest.realelems,
            ))
            copyto!(dest.data, realindices, bc.args[1].data, realindices)
        else
            copyto!(dest.data, bc.args[1])
        end
    else
        copyto!(dest.realdata, transform_broadcasted(bc, dest.data))
    end
    dest
end

"""
    begin_ghost_exchange!(Q::MPIStateArray; dependencies = nothing)

Begin the MPI halo exchange of the data stored in `Q`.  A KernelAbstractions
`Event` is returned that can be used as a dependency to end the exchange.
"""
function begin_ghost_exchange!(Q::MPIStateArray; dependencies = nothing)
    if !all(r -> r == MPI.REQUEST_NULL, Q.recvreq)
        error("The currrent ghost exchange must end before another begins.")
    end

    __Irecv!(Q)

    # wait on (prior) MPI sends
    @tic mpi_sendwait
    MPI.Waitall!(Q.sendreq)
    fill!(Q.sendreq, MPI.REQUEST_NULL)
    @toc mpi_sendwait

    @tic mpi_sendcopy
    stage = get_stage(Q.send_buffer)
    progress = () -> iprobe_and_yield(Q.mpicomm)
    event = fillsendbuf!(
        stage,
        Q.data,
        Q.vmapsend;
        dependencies = dependencies,
        progress = progress,
    )
    event = prepare_transfer!(
        Q.send_buffer;
        dependencies = event,
        progress = progress,
    )
    @toc mpi_sendcopy

    return event
end

"""
    end_ghost_exchange!(Q::MPIStateArray; dependencies = nothing)

This function blocks on the host until the ghost halo is received from MPI.  A
KernelAbstractions `Event` is returned that can be waited on to indicate when
the data is ready on the device.
"""
function end_ghost_exchange!(Q::MPIStateArray; dependencies = nothing)
    if any(r -> r == MPI.REQUEST_NULL, Q.recvreq)
        error("A ghost exchange must begin before it ends.")
    end

    progress = () -> iprobe_and_yield(Q.mpicomm)
    wait(CPU(), MultiEvent(dependencies), progress)

    __Isend!(Q)

    @tic mpi_recvwait
    MPI.Waitall!(Q.recvreq)
    fill!(Q.recvreq, MPI.REQUEST_NULL)
    @toc mpi_recvwait

    @tic mpi_recvcopy
    event = prepare_stage!(Q.recv_buffer; progress = progress)
    stage = get_stage(Q.recv_buffer)
    event = transferrecvbuf!(
        Q.data,
        stage,
        Q.vmaprecv;
        dependencies = event,
        progress = progress,
    )
    @toc mpi_recvcopy

    return event
end

function __Irecv!(Q)
    nnabr = length(Q.nabrtorank)
    transfer = get_transfer(Q.recv_buffer)

    for n in 1:nnabr
        # If this fails we haven't waited on previous recv!
        @assert Q.recvreq[n].buffer == nothing

        Q.recvreq[n] = MPI.Irecv!(
            (@view transfer[:, Q.nabrtovmaprecv[n]]),
            Q.nabrtorank[n],
            666,
            Q.mpicomm,
        )
    end
end

function __Isend!(Q)
    nnabr = length(Q.nabrtorank)
    transfer = get_transfer(Q.send_buffer)

    for n in 1:nnabr
        Q.sendreq[n] = MPI.Isend(
            (@view transfer[:, Q.nabrtovmapsend[n]]),
            Q.nabrtorank[n],
            666,
            Q.mpicomm,
        )
    end
end

function iprobe_and_yield(comm)
    MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, comm)
    yield()
end

# {{{ MPI Buffer handling
function fillsendbuf!(
    sendbuf,
    buf,
    vmapsend;
    dependencies = nothing,
    progress = yield,
)
    if length(vmapsend) == 0
        return MultiEvent(dependencies)
    end

    Np = size(buf, 1)
    nvar = size(buf, 2)

    event = kernel_fillsendbuf!(device(buf), 256)(
        Val(Np),
        Val(nvar),
        sendbuf,
        buf,
        vmapsend,
        length(vmapsend);
        ndrange = length(vmapsend),
        dependencies = dependencies,
        progress = progress,
    )

    return event
end

function transferrecvbuf!(
    buf,
    recvbuf,
    vmaprecv;
    dependencies = nothing,
    progress = yield,
)
    if length(vmaprecv) == 0
        return MultiEvent(dependencies)
    end

    Np = size(buf, 1)
    nvar = size(buf, 2)

    event = kernel_transferrecvbuf!(device(buf), 256)(
        Val(Np),
        Val(nvar),
        buf,
        recvbuf,
        vmaprecv,
        length(vmaprecv);
        ndrange = length(vmaprecv),
        dependencies = dependencies,
        progress = progress,
    )

    return event
end

# }}}

# Integral based metrics
function LinearAlgebra.norm(
    Q::MPIStateArray,
    p::Real = 2,
    weighted::Bool = true;
    dims = :,
)
    if weighted && ~isempty(Q.weights) && isfinite(p)
        W = @view Q.weights[:, :, Q.realelems]
        locnorm = weighted_norm_impl(Q.realdata, W, Val(p), dims)
    else
        locnorm = norm_impl(Q.realdata, Val(p), dims)
    end

    mpiop = isfinite(p) ? MPI.SUM : MPI.MAX
    if locnorm isa AbstractArray
        locnorm = convert(Array, locnorm)
    end
    @tic mpi_norm
    r = MPI.Allreduce(cpuify(locnorm), mpiop, Q.mpicomm)
    @toc mpi_norm
    isfinite(p) ? r .^ (1 // p) : r
end
LinearAlgebra.norm(Q::MPIStateArray, weighted::Bool; dims = :) =
    norm(Q, 2, weighted; dims = dims)

function LinearAlgebra.dot(
    Q1::MPIStateArray,
    Q2::MPIStateArray,
    weighted::Bool = true,
)
    @assert length(Q1.realdata) == length(Q2.realdata)

    if weighted && ~isempty(Q1.weights)
        W = @view Q1.weights[:, :, Q1.realelems]
        locnorm = weighted_dot_impl(Q1.realdata, Q2.realdata, W)
    else
        locnorm = dot_impl(Q1.realdata, Q2.realdata)
    end

    @tic mpi_dot
    r = MPI.Allreduce(locnorm, MPI.SUM, Q1.mpicomm)
    @toc mpi_dot
    return r
end

function euclidean_distance(A::MPIStateArray, B::MPIStateArray)
    # work around https://github.com/JuliaArrays/LazyArrays.jl/issues/66
    ArealQ = A.realdata
    BrealQ = B.realdata
    E = @~ (ArealQ .- BrealQ) .^ 2

    if ~isempty(A.weights)
        w = @view A.weights[:, :, A.realelems]
        E = @~ E .* w
    end

    locnorm = mapreduce(identity, +, E, init = zero(eltype(A)))
    @tic mpi_euclidean_distance
    r = sqrt(MPI.Allreduce(locnorm, MPI.SUM, A.mpicomm))
    @toc mpi_euclidean_distance
    return r
end

"""
    weightedsum(A[, states])

Compute the weighted sum of the `MPIStateArray` `A`. If `states` is specified on
the listed states are summed, otherwise all the states in `A` are used.

A typical use case for this is when the weights have been initialized with
quadrature weights from a grid, thus this becomes an integral approximation.
"""
function weightedsum(A::MPIStateArray, states = 1:size(A, 2))
    isempty(A.weights) && error("`weightedsum` requires weights")

    FT = eltype(A)
    states = SVector{length(states)}(states)

    C = @view A.data[:, states, A.realelems]
    w = @view A.weights[:, :, A.realelems]
    init = zero(DoubleFloat{FT})

    E = @~ DoubleFloat{FT}.(C) .* DoubleFloat{FT}.(w)

    locwsum = mapreduce(identity, +, E, init = init)

    @tic mpi_weightedsum
    # Need to use anomous function version of sum else MPI.jl using MPI_SUM
    r = FT(MPI.Allreduce(locwsum, (x, y) -> x + y, A.mpicomm))
    @toc mpi_weightedsum
    return r
end

# fast CPU local norm & dot implementations
function norm_impl(
    Q::SubArray{FT, N, A},
    ::Val{p},
    dims::Colon,
) where {FT, N, A <: Array, p}
    accum = isfinite(p) ? -zero(FT) : typemin(FT)
    @inbounds @simd for i in eachindex(Q)
        if isfinite(p)
            accum += abs(Q[i])^p
        else
            aQ_i = abs(Q[i])
            accum = ifelse(aQ_i > accum, aQ_i, accum)
        end
    end
    accum
end

function weighted_norm_impl(
    Q::SubArray{FT, N, A},
    W,
    ::Val{p},
    dims::Colon,
) where {FT, N, A <: Array, p}
    @assert isfinite(p)
    nq, ns, ne = size(Q)
    accum = -zero(FT)
    @inbounds for k in 1:ne, j in 1:ns
        @simd for i in 1:nq
            accum += W[i, 1, k] * abs(Q[i, j, k])^p
        end
    end
    accum
end

function dot_impl(
    Q1::SubArray{FT, N, A},
    Q2::SubArray{FT, N, A},
) where {FT, N, A <: Array}
    accum = -zero(FT)
    @inbounds @simd for i in eachindex(Q1)
        accum += Q1[i] * Q2[i]
    end
    accum
end

function weighted_dot_impl(
    Q1::SubArray{FT, N, A},
    Q2::SubArray{FT, N, A},
    W,
) where {FT, N, A <: Array}
    nq, ns, ne = size(Q1)
    accum = -zero(FT)
    @inbounds for k in 1:ne, j in 1:ns
        @simd for i in 1:nq
            accum += W[i, 1, k] * Q1[i, j, k] * Q2[i, j, k]
        end
    end
    accum
end

# GPU/generic local norm & dot implementations
function norm_impl(Q, ::Val{p}, dims = :) where {p}
    FT = eltype(Q)
    if !isfinite(p)
        f, op = abs, max
    elseif p == 1
        f, op = abs, +
    elseif p == 2
        f, op = abs2, +
    else
        f, op = x -> abs(x)^p, +
    end
    mapreduce(f, op, Q, dims = dims)
end

function weighted_norm_impl(Q, W, ::Val{p}, dims = :) where {p}
    @assert isfinite(p)
    FT = eltype(Q)
    if p == 1
        E = @~ @. W * abs(Q)
    elseif p == 2
        E = @~ @. W * abs2(Q)
    else
        E = @~ @. W * abs(Q)^p
    end
    op, init = +, zero(FT)
    reduce(op, E, init = init, dims = dims)
end

function dot_impl(Q1, Q2)
    FT = eltype(Q1)
    E = @~ @. Q1 * Q2
    mapreduce(identity, +, E, init = zero(FT))
end

weighted_dot_impl(Q1, Q2, W) = dot_impl(@~ @. W * Q1, Q2)

function Base.mapreduce(f, op, Q::MPIStateArray; kw...)
    locreduce = mapreduce(f, op, realview(Q); kw...)
    MPI.Allreduce(cpuify(locreduce), op, Q.mpicomm)
end

# Arrays and CuArrays have different reduction machinery
# until we can figure this out, add special cases to make common functions work
function Base.mapreduce(
    ::typeof(identity),
    ::Union{typeof(+), typeof(Base.add_sum)},
    Q::MPIStateArray;
    kw...,
)
    locreduce = sum(realview(Q); kw...)
    MPI.Allreduce(cpuify(locreduce), +, Q.mpicomm)
end
function Base.mapreduce(
    ::typeof(identity),
    ::typeof(min),
    Q::MPIStateArray;
    kw...,
)
    locreduce = minimum(realview(Q); kw...)
    MPI.Allreduce(cpuify(locreduce), min, Q.mpicomm)
end
function Base.mapreduce(
    ::typeof(identity),
    ::typeof(max),
    Q::MPIStateArray;
    kw...,
)
    locreduce = maximum(realview(Q); kw...)
    MPI.Allreduce(cpuify(locreduce), max, Q.mpicomm)
end

# `realview` and `device` are helpers that enable
# testing ODESolvers and LinearSolvers without using MPIStateArrays
# They could be potentially useful elsewhere and exported but probably need
# better names, for example `device` is also defined in CUDAdrv

device(::Union{Array, SArray, MArray}) = CPU()
device(Q::MPIStateArray) = device(Q.data)

realview(Q::Union{Array, SArray, MArray}) = Q
realview(Q::MPIStateArray) = Q.realdata

device(::CuArray) = CUDA()
realview(Q::CuArray) = Q

# transform all arguments of `bc` from MPIStateArrays to CuArrays
# and replace CPU function with GPU variants
function transform_broadcasted(bc::Broadcasted, ::CuArray)
    transform_cuarray(bc)
end
function transform_cuarray(bc::Broadcasted)
    Broadcasted(CuArrays.cufunc(bc.f), transform_cuarray.(bc.args), bc.axes)
end
transform_cuarray(mpisa::MPIStateArray) = mpisa.realdata
transform_cuarray(x) = x

# @init tictoc()

using KernelAbstractions.Extras: @unroll

@kernel function kernel_fillsendbuf!(
    ::Val{Np},
    ::Val{nvar},
    sendbuf,
    buf,
    vmapsend,
    nvmapsend,
) where {Np, nvar}

    i = @index(Global, Linear)
    @inbounds begin
        e, n = fldmod1(vmapsend[i], Np)
        @unroll for s in 1:nvar
            sendbuf[s, i] = buf[n, s, e]
        end
    end
end

@kernel function kernel_transferrecvbuf!(
    ::Val{Np},
    ::Val{nvar},
    buf,
    recvbuf,
    vmaprecv,
    nvmaprecv,
) where {Np, nvar}

    i = @index(Global, Linear)
    @inbounds begin
        e, n = fldmod1(vmaprecv[i], Np)
        @unroll for s in 1:nvar
            buf[n, s, e] = recvbuf[s, i]
        end
    end
end

end


text:
function _permutedims(S::FourthOrderTensor{dim}, idx::NTuple{4,Int}) where dim
    sort([idx...]) == [1,2,3,4] || throw(ArgumentError("Missing index."))
    neworder = sortperm([idx...])
    f = (i,j,k,l) -> S[[i,j,k,l][neworder]...]
    return Tensor{4,dim}(f)
end

@testsection "tensor ops" begin
for T in (Float32, Float64, F64), dim in (1,2,3)
println("T = $T, dim = $dim")
AA = rand(Tensor{4, dim, T})
BB = rand(Tensor{4, dim, T})
A = rand(Tensor{2, dim, T})
B = rand(Tensor{2, dim, T})
a = rand(Tensor{1, dim, T})
b = rand(Tensor{1, dim, T})

AA_sym = rand(SymmetricTensor{4, dim, T})
BB_sym = rand(SymmetricTensor{4, dim, T})
A_sym = rand(SymmetricTensor{2, dim, T})
B_sym = rand(SymmetricTensor{2, dim, T})

i,j,k,l = rand(1:dim,4)

@testsection "double contraction" begin
    # 4 - 4
    @test vec((@inferred dcontract(AA, BB))::Tensor{4, dim, T})                  ≈ vec(collect(reshape(vec(AA), (dim^2, dim^2))) * collect(reshape(vec(BB), (dim^2, dim^2))))
    @test vec((@inferred dcontract(AA_sym, BB))::Tensor{4, dim, T})              ≈ vec(collect(reshape(vec(AA_sym), (dim^2, dim^2))) * collect(reshape(vec(BB), (dim^2, dim^2))))
    @test vec((@inferred dcontract(AA, BB_sym))::Tensor{4, dim, T})              ≈ vec(collect(reshape(vec(AA), (dim^2, dim^2))) * collect(reshape(vec(BB_sym), (dim^2, dim^2))))
    @test vec((@inferred dcontract(AA_sym, BB_sym))::SymmetricTensor{4, dim, T}) ≈ vec(collect(reshape(vec(AA_sym), (dim^2, dim^2))) * collect(reshape(vec(BB_sym), (dim^2, dim^2))))
    @test dcontract(convert(Tensor, AA_sym), convert(Tensor, BB_sym))            ≈ dcontract(AA_sym, BB_sym)

    # 2 - 4
    @test (@inferred dcontract(AA, A))::Tensor{2, dim, T}                  ≈ reshape(collect(reshape(vec(AA), (dim^2, dim^2))) * collect(reshape(vec(A), (dim^2,))), dim, dim)
    @test (@inferred dcontract(AA_sym, A))::SymmetricTensor{2, dim, T}     ≈ reshape(collect(reshape(vec(AA_sym), (dim^2, dim^2))) * collect(reshape(vec(A), (dim^2,))), dim, dim)
    @test (@inferred dcontract(AA, A_sym))::Tensor{2, dim, T}              ≈ reshape(collect(reshape(vec(AA), (dim^2, dim^2))) * collect(reshape(vec(A_sym), (dim^2,))), dim, dim)
    @test (@inferred dcontract(AA_sym, A_sym))::SymmetricTensor{2, dim, T} ≈ reshape(collect(reshape(vec(AA_sym), (dim^2, dim^2))) * collect(reshape(vec(A_sym), (dim^2,))), dim, dim)
    @test (@inferred dcontract(A, AA))::Tensor{2, dim, T}                  ≈ reshape(collect(reshape(vec(AA), (dim^2, dim^2))') * collect(reshape(vec(A), (dim^2,))), dim, dim)
    @test (@inferred dcontract(A_sym, AA))::Tensor{2, dim, T}              ≈ reshape(collect(reshape(vec(AA), (dim^2, dim^2))') * collect(reshape(vec(A_sym), (dim^2,))), dim, dim)
    @test (@inferred dcontract(A, AA_sym))::SymmetricTensor{2, dim, T}     ≈ reshape(collect(reshape(vec(AA_sym), (dim^2, dim^2))') * collect(reshape(vec(A), (dim^2,))), dim, dim)
    @test (@inferred dcontract(A_sym, AA_sym))::SymmetricTensor{2, dim, T} ≈ reshape(collect(reshape(vec(AA_sym), (dim^2, dim^2))') * collect(reshape(vec(A_sym), (dim^2,))), dim, dim)
    @test dcontract(convert(Tensor, AA_sym), convert(Tensor, A_sym))       ≈ dcontract(AA_sym, A_sym)

    # 2 - 2
    @test (@inferred dcontract(A, B))::T         ≈ sum(vec(A) .* vec(B))
    @test (@inferred dcontract(A_sym, B))::T     ≈ sum(vec(A_sym) .* vec(B))
    @test (@inferred dcontract(A, B_sym))::T     ≈ sum(vec(A) .* vec(B_sym))
    @test (@inferred dcontract(A_sym, B_sym))::T ≈ sum(vec(A_sym) .* vec(B_sym))
end # of testsection

@testsection "outer products (otimes, otimesu, otimesl)" begin
    # binary otimes
    @test             (@inferred otimes(a, b))::Tensor{2, dim, T}                                  ≈ Array(a) * Array(b)'
    @test reshape(vec((@inferred otimes(A, B))::Tensor{4, dim, T}), dim^2, dim^2)                  ≈ vec(A) * vec(B)'
    @test reshape(vec((@inferred otimes(A_sym, B))::Tensor{4, dim, T}), dim^2, dim^2)              ≈ vec(A_sym) * vec(B)'
    @test reshape(vec((@inferred otimes(A, B_sym))::Tensor{4, dim, T}), dim^2, dim^2)              ≈ vec(A) * vec(B_sym)'
    @test reshape(vec((@inferred otimes(A_sym, B_sym))::SymmetricTensor{4, dim, T}), dim^2, dim^2) ≈ vec(A_sym) * vec(B_sym)'
    # unary otimes
    @test (@inferred otimes(a))::SymmetricTensor{2, dim, T} ≈ otimes(a, a)

    # binary otimesu
    @test (@inferred otimesu(A, B))::Tensor{4, dim, T}         ≈ _permutedims(otimes(A, B), (1,3,2,4))
    @test (@inferred otimesu(A_sym, B))::Tensor{4, dim, T}     ≈ _permutedims(otimes(A_sym, B), (1,3,2,4))
    @test (@inferred otimesu(A, B_sym))::Tensor{4, dim, T}     ≈ _permutedims(otimes(A, B_sym), (1,3,2,4))
    @test (@inferred otimesu(A_sym, B_sym))::Tensor{4, dim, T} ≈ _permutedims(otimes(A_sym, B_sym), (1,3,2,4))

    # binary otimesl
    @test (@inferred otimesl(A, B))::Tensor{4, dim, T}         ≈ _permutedims(otimes(A, B), (1,3,4,2))
    @test (@inferred otimesl(A_sym, B))::Tensor{4, dim, T}     ≈ _permutedims(otimes(A_sym, B), (1,3,4,2))
    @test (@inferred otimesl(A, B_sym))::Tensor{4, dim, T}     ≈ _permutedims(otimes(A, B_sym), (1,3,4,2))
    @test (@inferred otimesl(A_sym, B_sym))::Tensor{4, dim, T} ≈ _permutedims(otimes(A_sym, B_sym), (1,3,4,2))
end # of testsection

@testsection "dot products" begin
    # 1 - 1
    @test (@inferred dot(a, b))::T               ≈ sum(Array(a) .* Array(b))
    # 1 - 2
    @test (@inferred dot(A, b))::Vec{dim, T}     ≈ Array(A) * Array(b)
    @test (@inferred dot(A_sym, b))::Vec{dim, T} ≈ Array(A_sym) * Array(b)
    @test (@inferred dot(a, B))::Vec{dim, T}     ≈ Array(B)' * Array(a)
    @test (@inferred dot(a, B_sym))::Vec{dim, T} ≈ Array(B_sym)' * Array(a)
    # 2 - 2
    # binary
    @test (@inferred dot(A, B))::Tensor{2, dim, T}         ≈ Array(A) * Array(B)
    @test (@inferred dot(A_sym, B))::Tensor{2, dim, T}     ≈ Array(A_sym) * Array(B)
    @test (@inferred dot(A, B_sym))::Tensor{2, dim, T}     ≈ Array(A) * Array(B_sym)
    @test (@inferred dot(A_sym, B_sym))::Tensor{2, dim, T} ≈ Array(A_sym) * Array(B_sym)
    # unary
    @test (@inferred dot(A_sym))::SymmetricTensor{2, dim, T}  ≈ dot(A_sym, A_sym)
    @test (@inferred tdot(A))::SymmetricTensor{2, dim, T}     ≈ dot(transpose(A), A)
    @test (@inferred tdot(A_sym))::SymmetricTensor{2, dim, T} ≈ dot(transpose(A_sym), A_sym)
    @test (@inferred dott(A))::SymmetricTensor{2, dim, T}     ≈ dot(A, transpose(A))
    @test (@inferred dott(A_sym))::SymmetricTensor{2, dim, T} ≈ dot(A_sym, transpose(A_sym))
    # 2 - 4
    @test (@inferred dot(AA, B))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(AA), (dim^3, dim))) * collect(reshape(vec(B), (dim, dim))), (dim, dim, dim, dim))
    @test (@inferred dot(B, AA))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(B), (dim, dim))) * collect(reshape(vec(AA), (dim, dim^3))), (dim, dim, dim, dim))
    @test (@inferred dot(AA_sym, B))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(AA_sym), (dim^3, dim))) * collect(reshape(vec(B), (dim, dim))), (dim, dim, dim, dim))
    @test (@inferred dot(B, AA_sym))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(B), (dim, dim))) * collect(reshape(vec(AA_sym), (dim, dim^3))), (dim, dim, dim, dim))
    @test (@inferred dot(AA, B_sym))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(AA), (dim^3, dim))) * collect(reshape(vec(B_sym), (dim, dim))), (dim, dim, dim, dim))
    @test (@inferred dot(B_sym, AA))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(B_sym), (dim, dim))) * collect(reshape(vec(AA), (dim, dim^3))), (dim, dim, dim, dim))
    @test (@inferred dot(AA_sym, B_sym))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(AA_sym), (dim^3, dim))) * collect(reshape(vec(B_sym), (dim, dim))), (dim, dim, dim, dim))
    @test (@inferred dot(B_sym, AA_sym))::Tensor{4, dim, T} ≈ reshape(collect(reshape(vec(B_sym), (dim, dim))) * collect(reshape(vec(AA_sym), (dim, dim^3))), (dim, dim, dim, dim))
end # of testsection

@testsection "symmetric/skew-symmetric" begin
    if dim == 1 # non-symmetric tensors are symmetric
        @test (@inferred issymmetric(A))
        @test (@inferred issymmetric(AA))
    elseif dim != 1
        @test !(@inferred issymmetric(A))
        @test !(@inferred issymmetric(AA))
        @test !(@inferred ismajorsymmetric(AA))
        @test !(@inferred isminorsymmetric(AA))
        @test_throws InexactError convert(typeof(A_sym),A)
        @test_throws InexactError convert(typeof(AA_sym),AA)
    end
    @test (@inferred issymmetric(A_sym))
    @test (@inferred issymmetric(AA_sym))
    @test (@inferred isminorsymmetric(AA_sym))
    @test (@inferred issymmetric(symmetric(A)))
    @test (@inferred issymmetric(A + A'))

    @test (@inferred symmetric(A))::SymmetricTensor{2, dim, T}     ≈ 0.5(A + A')
    @test (@inferred symmetric(A_sym))::SymmetricTensor{2, dim, T} ≈ A_sym
    @test convert(typeof(A_sym),convert(Tensor,symmetric(A))) ≈ symmetric(A)

    @test (@inferred symmetric(AA))::SymmetricTensor{4, dim, T} ≈ (@inferred minorsymmetric(AA))::SymmetricTensor{4, dim, T}
    @test minorsymmetric(AA)[i,j,k,l] ≈ minorsymmetric(AA)[j,i,l,k]
    @test issymmetric(convert(Tensor,minorsymmetric(AA)))
    @test isminorsymmetric(convert(Tensor,minorsymmetric(AA)))
    @test (@inferred symmetric(AA_sym))::SymmetricTensor{4, dim, T} ≈ (@inferred minorsymmetric(AA_sym))::SymmetricTensor{4, dim, T}
    @test minorsymmetric(AA_sym)[i,j,k,l] ≈ minorsymmetric(AA_sym)[j,i,l,k]
    @test minorsymmetric(AA_sym) ≈ AA_sym
    @test issymmetric(convert(Tensor,minorsymmetric(AA_sym)))
    @test isminorsymmetric(convert(Tensor,minorsymmetric(AA_sym)))

    @test convert(typeof(AA_sym),convert(Tensor,minorsymmetric(AA))) ≈ minorsymmetric(AA)

    @test ((@inferred majorsymmetric(AA))::Tensor{4, dim, T})[i,j,k,l] ≈ 0.5*(AA[i,j,k,l] + AA[k,l,i,j])
    @test majorsymmetric(AA)[i,j,k,l] ≈ majorsymmetric(AA)[k,l,i,j]
    @test ismajorsymmetric(majorsymmetric(AA))
    @test ((@inferred majorsymmetric(AA_sym))::SymmetricTensor{4, dim, T})[i,j,k,l] ≈ 0.5*(AA_sym[i,j,k,l] + AA_sym[k,l,i,j])
    @test majorsymmetric(AA_sym)[i,j,k,l] ≈ majorsymmetric(AA_sym)[k,l,i,j]
    @test ismajorsymmetric(majorsymmetric(AA_sym))

    @test (@inferred skew(A))::Tensor{2, dim, T} ≈ 0.5(A - A')
    @test (@inferred skew(A_sym))::Tensor{2, dim, T} ≈ zero(A_sym)

    # Identities
    @test A ≈ symmetric(A) + skew(A)
    @test skew(A) ≈ -skew(A)'
    @test tr(skew(A)) ≈ 0.0
    @test tr(symmetric(A)) ≈ tr(A)
end # of testsection

@testsection "transpose" begin
    @test (@inferred transpose(A))::Tensor{2, dim, T} ≈ Array(A)'
    @test transpose(transpose(A)) ≈ A
    @test (@inferred transpose(A_sym))::SymmetricTensor{2, dim, T} ≈ A_sym ≈ Array(A_sym)'
    @test transpose(transpose(A_sym)) ≈ A_sym

    @test (@inferred transpose(AA))::Tensor{4, dim, T} ≈ (@inferred minortranspose(AA))::Tensor{4, dim, T}
    @test AA[i,j,k,l] ≈ minortranspose(AA)[j,i,l,k]
    @test AA_sym[i,j,k,l] ≈ minortranspose(AA_sym)[j,i,l,k]
    @test (@inferred transpose(AA_sym))::SymmetricTensor{4, dim, T} ≈ (@inferred minortranspose(AA_sym))::SymmetricTensor{4, dim, T} ≈ AA_sym
    @test minortranspose(minortranspose(AA)) ≈ AA
    @test minortranspose(minortranspose(AA_sym)) ≈ AA_sym

    @test majortranspose((@inferred majortranspose(AA))::Tensor{4, dim, T}) ≈ AA
    @test majortranspose((@inferred majortranspose(AA_sym))::Tensor{4, dim, T}) ≈ AA_sym
    @test AA[i,j,k,l] ≈ majortranspose(AA)[k,l,i,j]
    @test AA_sym[i,j,k,l] ≈ majortranspose(AA_sym)[k,l,i,j]

    @test minortranspose(AA) ≈ _permutedims(AA,(2,1,4,3))
    @test minortranspose(AA_sym) ≈ _permutedims(AA_sym,(2,1,4,3))
    @test majortranspose(AA) ≈ _permutedims(AA,(3,4,1,2))
    @test majortranspose(AA_sym) ≈ _permutedims(AA_sym,(3,4,1,2))
end # of testsection

@testsection "cross product" begin
    @test (@inferred a × a)::Vec{3, T} ≈ Vec{3, T}((0.0,0.0,0.0))
    @test a × b ≈ -b × a
    if dim == 2
        ad = Vec{2, T}((1.0,0.0))
        ad2 = Vec{2, T}((0.0,1.0))
        @test (@inferred ad × ad2)::Vec{3, T} ≈ Vec{3, T}((0.0, 0.0, 1.0))
    end
    if dim == 3
        ad = Vec{3, T}((1.0,0.0,0.0))
        ad2 = Vec{3, T}((0.0,1.0,0.0))
        @test (@inferred ad × ad2)::Vec{3, T} ≈ Vec{3, T}((0.0, 0.0, 1.0))
    end
    if T == Float64 # mixed eltype
        @test rand(Vec{dim,Float64}) × rand(Vec{dim,Float32}) isa Vec{3,Float64}
    end
end # of testsection

@testsection "special" begin
    AAT = Tensor{4, dim, T}((i,j,k,l) -> AA_sym[i,l,k,j])
    @test AAT ⊡ (b ⊗ a) ≈ (@inferred dotdot(a, AA_sym, b))::Tensor{2, dim, T}
end # of testsection

if dim == 3
    @testsection "rotation" begin
        x = eᵢ(Vec{3, T}, 1)
        y = eᵢ(Vec{3, T}, 2)
        z = eᵢ(Vec{3, T}, 3)

        @test (@inferred rotate(z, z, rand(T)))::Vec{3, T} ≈ z
        @test rotate(2*z, y, π/2) ≈ 2*x
        @test rotate(3*z, y, π) ≈ -3*z
        @test rotate(x+y+z, z, π/4) ≈ Vec{3}((0.0,√2,1.0))

        @test rotate(a, b, 0) ≈ a
        @test rotate(a, b, π) ≈ rotate(a, b, -π)
        @test rotate(a, b, π/2) ≈ rotate(a, -b, -π/2)

        @test rotate(A, x, π) ≈ rotate(A, x, -π)
        @test rotate(rotate(rotate(A, x, π), y, π), z, π) ≈ A
        @test rotate(A, a, 0) ≈ A
        @test rotate(A, a, π/2) ≈ rotate(A, -a, -π/2)
        A_sym_T = convert(Tensor{2,3}, A_sym)
        @test rotate(A_sym, x, π)::SymmetricTensor ≈ rotate(A_sym, x, -π)::SymmetricTensor ≈ rotate(A_sym_T, x, -π)
        @test rotate(rotate(rotate(A_sym, x, π), y, π), z, π)::SymmetricTensor ≈ A_sym
        @test rotate(A_sym, a, 0) ≈ A_sym
        @test rotate(A_sym, a, π/2) ≈ rotate(A_sym, -a, -π/2)

        @test rotate(AA, x, π)::Tensor ≈ rotate(AA, x, -π)
        @test rotate(rotate(rotate(AA, x, π), y, π), z, π) ≈ AA
        @test rotate(AA, a, 0) ≈ AA
        @test rotate(AA, a, π/2) ≈ rotate(AA, -a, -π/2)
        AA_sym_T = convert(Tensor{4,3}, AA_sym)
        @test rotate(AA_sym, x, π)::SymmetricTensor ≈ rotate(AA_sym, x, -π) ≈ rotate(AA_sym_T, x, π)
        @test rotate(rotate(rotate(AA_sym, x, π), y, π), z, π) ≈ AA_sym
        @test rotate(AA_sym, a, 0) ≈ AA_sym
        @test rotate(AA_sym, a, π/2) ≈ rotate(AA_sym, -a, -π/2)

        v1, v2, v3, v4, axis = [rand(Vec{3,T}) for _ in 1:5]
        α = rand(T) * π
        R = Tensors.rotation_matrix(axis, α)
        v1v2v3v4 = (v1 ⊗ v2) ⊗ (v3 ⊗ v4)
        Rv1v2v3v4 = ((R ⋅ v1) ⊗ (R ⋅ v2)) ⊗ ((R ⋅ v3) ⊗ (R ⋅ v4))
        @test rotate(v1v2v3v4, axis, α) ≈ Rv1v2v3v4
        v1v1v2v2 = otimes(v1) ⊗ otimes(v2)
        Rv1v1v2v2 = otimes(R ⋅ v1) ⊗ otimes(R ⋅ v2)
        @test rotate(v1v1v2v2, axis, α) ≈ Rv1v1v2v2
    end
end

@testsection "tovoigt/fromvoigt" begin
    # https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/books/usb/default.htm?startat=pt01ch01s02aus02.html
    abaqus = dim == 1 ? fill(1, 1, 1) : dim == 2 ? [1 3; 0 2] : [1 4 6; 0 2 5; 0 0 3]
    T2 = T(2)
    @test (@inferred tovoigt(AA)) * (@inferred tovoigt(A)) ≈ tovoigt(AA ⊡ A)
    @test (@inferred tovoigt(AA_sym)) * (@inferred tovoigt(A_sym, offdiagscale=T2)) ≈ tovoigt(AA_sym ⊡ A_sym)

    nc(x::SecondOrderTensor) = length(x.data)
    nc(x::FourthOrderTensor) = Int(sqrt(length(x.data)))
    x = zeros(T, nc(A) + 1); tovoigt!(x, A; offset=1)
    @test x[2:end] == tovoigt(A)
    x = zeros(T, nc(A_sym) + 1); tovoigt!(x, A_sym; offset=1)
    @test x[2:end] == tovoigt(A_sym)
    x = zeros(T, nc(AA) + 1, nc(AA) + 1); tovoigt!(x, AA; offset_i=1, offset_j=1)
    @test x[2:end, 2:end] == tovoigt(AA)
    x = zeros(T, nc(AA_sym) + 1, nc(AA_sym) + 1); tovoigt!(x, AA_sym; offset_i=1, offset_j=1)
    @test x[2:end, 2:end] == tovoigt(AA_sym)

    x = zeros(T, nc(A) + 1); tomandel!(x, A; offset=1)
    @test x[2:end] == tomandel(A)
    x = zeros(T, nc(A_sym) + 1); tomandel!(x, A_sym; offset=1)
    @test x[2:end] == tomandel(A_sym)
    x = zeros(T, nc(AA) + 1, nc(AA) + 1); tomandel!(x, AA; offset_i=1, offset_j=1)
    @test x[2:end, 2:end] == tomandel(AA)
    x = zeros(T, nc(AA_sym) + 1, nc(AA_sym) + 1); tomandel!(x, AA_sym; offset_i=1, offset_j=1)
    @test x[2:end, 2:end] == tomandel(AA_sym)

    @test (@inferred fromvoigt(Tensor{2,dim}, tovoigt(A))) ≈ A
    @test (@inferred fromvoigt(Tensor{4,dim}, tovoigt(AA))) ≈ AA
    @test (@inferred fromvoigt(SymmetricTensor{2,dim}, tovoigt(A_sym, offdiagscale=T2), offdiagscale=T2)) ≈ A_sym
    @test (@inferred fromvoigt(SymmetricTensor{4,dim}, tovoigt(AA_sym, offdiagscale=T2), offdiagscale=T2)) ≈ AA_sym

    @test (@inferred fromvoigt(SymmetricTensor{2,dim}, tovoigt(A_sym; order=abaqus); order=abaqus)) ≈ A_sym
    @test (@inferred fromvoigt(SymmetricTensor{4,dim}, tovoigt(AA_sym; order=abaqus); order=abaqus)) ≈ AA_sym

    @test (@inferred tomandel(AA_sym)) * (@inferred tomandel(A_sym)) ≈ tomandel(AA_sym ⊡ A_sym)
    @test (@inferred frommandel(SymmetricTensor{2,dim}, tomandel(A_sym))) ≈ A_sym
    @test (@inferred frommandel(SymmetricTensor{4,dim}, tomandel(AA_sym))) ≈ AA_sym
    @test (@inferred tomandel(AA)) * (@inferred tomandel(A)) ≈ tomandel(AA ⊡ A)
    @test (@inferred frommandel(Tensor{2,dim}, tomandel(A))) ≈ A
    @test (@inferred frommandel(Tensor{4,dim}, tomandel(AA))) ≈ AA

    if T==Float64
        num_components = Int((dim^2+dim)/2)
        @test isa(fromvoigt(SymmetricTensor{2,dim}, rand(num_components) .> 0.5), SymmetricTensor{2,dim,Bool})
        @test isa(fromvoigt(SymmetricTensor{4,dim}, rand(num_components,num_components) .> 0.5), SymmetricTensor{4,dim,Bool})
    end
end
end # of testsection
end # of testsection


text:
#First delete any existing nc-files
fn1=tempname()
fn2=tempname()
fn3=tempname()

# Test Medium level Interface
# Test Dimension Creation
d1 = NcDim("Dim1",2;values=[5.0,10.0],atts=@Compat.AnyDict("units"=>"deg C"));
d2 = NcDim("Dim2",collect(1:10));
d3 = NcDim("Dim3",20;atts=@Compat.AnyDict("max"=>10));

# Test Variable creation
v1 = NcVar("v1",[d1,d2,d3],compress=5) 						# With several dims in an Array, and compressed
v2 = NcVar("v2",[d1,d2,d3],atts=@Compat.AnyDict("a1"=>"varatts"))  # with given attributes
v3 = NcVar("v3",d1)
vs = NcVar("vstr",d2,t=String) 								# with a single dimension
tlist = [Float64, Float32, Int32, Int16, Int8]
vt = Array(NcVar, length(tlist))
for i= 1:length(tlist)
  vt[i]=NcVar("vt$i",d2,t = tlist[i])
end

# Creating Files
nc1 = NetCDF.create(fn1,v1,vs,mode=NC_NETCDF4);
nc2 = NetCDF.create(fn2,v2,v3,gatts=@Compat.AnyDict("Some global attributes"=>2010),mode=NC_64BIT_OFFSET);
nc3 = NetCDF.create(fn3,vt,mode=NC_CLASSIC_MODEL);

#Test Adding attributes
NetCDF.putatt(nc1,"v1",@Compat.Dict("Additional String attribute"=>"att"))
NetCDF.putatt(nc1,"global",@Compat.Dict("Additional global attribute"=>"gatt"))
NetCDF.putatt(nc1,"global",@Compat.Dict("Additional Int8 attribute"=>Int8(20),
							"Additional Int16 attribute"=>Int16(20),
							"Additional Int32 attribute"=>Int32(20),
							"Additional Float32 attribute"=>Float32(20),
							"Additional Float64 attribute"=>Float64(20)))
NetCDF.putatt(nc1,"v1",    @Compat.Dict("Additional Int8 array attribute"=>Int8[i for i in 1:20],
							"Additional Int16 array attribute"=>Int16[i for i in 1:20],
							"Additional Int32 array attribute"=>Int32[i for i in 1:20],
							"Additional Float32 array attribute"=>Float32[i for i in 1:20],
							"Additional Float64 array attribute"=>Float64[i for i in 1:20]))

#Test writing data
#
#First generate the data
x1 = rand(2,10,20)
x2 = rand(2,10,20)
xt=Array(Any,length(tlist))
for i=1:length(tlist)
    xt[i]=rand(tlist[i],10)
end
xs=Array(AbstractString,10)
xs[1]="a"
xs[2]="bb"
xs[3]="ccc"
xs[4]="dddd"
xs[5]="eeeee"
xs[6]="ffffff"
xs[7]="ggggggg"
xs[8]="hhhhhhhh"
xs[9]="iiiiiiiii"
xs[10]="jjjjjjjjjj"
#
# And write it
#
NetCDF.putvar(nc1,"v1",x1)
NetCDF.putvar(nc1,"vstr",xs)
#Test sequential writing along one dimension
for i=1:10
  NetCDF.putvar(nc2,"v2",x2[:,i,:],start=[1,i,1],count=[-1,1,-1])
end
#Test automatic type conversion
x4=[1,2]
NetCDF.putvar(nc2,"v3",x4)

for i=1:length(tlist)
  NetCDF.putvar(nc3,"vt$i",xt[i])
end

NetCDF.close(nc1)
NetCDF.close(nc2)
NetCDF.close(nc3)

##Read the data back and compare
nc1 = NetCDF.open(fn1,mode=NC_NOWRITE);
nc2 = NetCDF.open(fn2,mode=NC_NOWRITE);
nc3 = NetCDF.open(fn3,mode=NC_NOWRITE);

@test x1==NetCDF.readvar(nc1,"v1")
@test xs==NetCDF.readvar(nc1,"vstr")
@test x2==NetCDF.readvar(nc2,"v2")
@test x4==NetCDF.readvar(nc2,"v3")


text:
function erdos_renyi_test()
    @test_throws DomainError erdos_renyi_undirected(10,11.)
    @test_throws DomainError erdos_renyi_directed(10,11.)

    n = 100
    avgdegs = linspace(1.,2*log(n),100)
    compsizes = map( (dbar) ->
            maximum(scomponents(erdos_renyi_undirected(n,dbar)).sizes),
        avgdegs )

    @test is_undirected(erdos_renyi_undirected(10,2.))

    erdos_renyi_undirected(5, 0.)
    erdos_renyi_undirected(5, 1.)
    @test all(diag(sparse_transpose(erdos_renyi_undirected(10, 0.5))) .== 0.)
    erdos_renyi_directed(5, 0.)
    erdos_renyi_directed(5, 1.)
    @test all(diag(sparse_transpose(erdos_renyi_directed(10, 0.5))) .== 0.)

end

function chung_lu_test()
    @test_throws ArgumentError chung_lu_undirected([1,2],2)
    @test_throws ArgumentError chung_lu_undirected([1,3])
    @test size(chung_lu_undirected([1,2]),1) == 2
end

function havel_hakimi_test()
    @test is_graphical_sequence([1,1]) == true
    @test is_graphical_sequence([0,1,1]) == true
    @test is_graphical_sequence([2,2,2]) == true
    @test is_graphical_sequence([1,1,2]) == true
    @test is_graphical_sequence([1,2,2]) == false
    @test is_graphical_sequence([1,1,3]) == false
    @test is_graphical_sequence([0]) == true
    @test is_graphical_sequence(Int[]) == true
    @test is_graphical_sequence([1,1,1,1]) == true

    @test full(sparse_transpose(havel_hakimi_graph(Int[]))) == zeros(Bool,0,0)
    @test full(sparse_transpose(havel_hakimi_graph([0]))) == zeros(Bool,1,1)
    @test full(sparse_transpose(havel_hakimi_graph([1,1]))) == [0 1; 1 0]
    @test full(sparse_transpose(havel_hakimi_graph([2,2,2]))) == [0 1 1; 1 0 1; 1 1 0]
    six_star = [0 1 1 1 1 1; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0]
    @test full(sparse_transpose(havel_hakimi_graph([5,1,1,1,1,1]))) == six_star

   @test_throws ArgumentError is_graphical_sequence([-1])
   @test_throws ArgumentError havel_hakimi_graph([1,2,2])

end

function generators_test()

erdos_renyi_test()
chung_lu_test()
havel_hakimi_test()
return true

end


text:
# ---
# title: 185. Department Top Three Salaries
# id: problem185
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/department-top-three-salaries/description/>
# hidden: true
# ---
#
# The `Employee` table holds all employees. Every employee has an Id, and there
# is also a column for the department Id.
#
#
#
#     +----+-------+--------+--------------+
#     | Id | Name  | Salary | DepartmentId |
#     +----+-------+--------+--------------+
#     | 1  | Joe   | 85000  | 1            |
#     | 2  | Henry | 80000  | 2            |
#     | 3  | Sam   | 60000  | 2            |
#     | 4  | Max   | 90000  | 1            |
#     | 5  | Janet | 69000  | 1            |
#     | 6  | Randy | 85000  | 1            |
#     | 7  | Will  | 70000  | 1            |
#     +----+-------+--------+--------------+
#
#
# The `Department` table holds all departments of the company.
#
#
#
#     +----+----------+
#     | Id | Name     |
#     +----+----------+
#     | 1  | IT       |
#     | 2  | Sales    |
#     +----+----------+
#
#
# Write a SQL query to find employees who earn the top three salaries in each of
# the department. For the above tables, your SQL query should return the
# following rows (order of rows does not matter).
#
#
#
#     +------------+----------+--------+
#     | Department | Employee | Salary |
#     +------------+----------+--------+
#     | IT         | Max      | 90000  |
#     | IT         | Randy    | 85000  |
#     | IT         | Joe      | 85000  |
#     | IT         | Will     | 70000  |
#     | Sales      | Henry    | 80000  |
#     | Sales      | Sam      | 60000  |
#     +------------+----------+--------+
#
#
# Explanation:
#
# In IT department, Max earns the highest salary, both Randy and Joe earn the
# second highest salary, and Will earns the third highest salary. There are only
# two employees in the Sales department, Henry earns the highest salary while
# Sam earns the second highest salary.
#
#
## @lc code=start
using LeetCode

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


text:
densities_T = [25.0, 38.0, 1.0] # TCR, pMHC, Ac
densities_TC1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_TC2 = [8.0,  38.0, 17.0, 1.0]
densities_TC3 = [3.0,  38.0, 22.0, 1.0]

rates_T = [-3.37638358011457
            -0.302013187183543]
rates_TC1 = [-1.72689171836867
            -0.646695212785203]
rates_TC2 = [-0.644293816490510
            0.617796010276307]
rates_TC3 = [0.891017383250648
            2.16934160248305]

sd_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/n_stds_for_julia.csv")
nT_sds = sd_data[:,1]
nTC1_sds = Array(sd_data[:,2:3])
nTC2_sds = Array(sd_data[:,4:5])
nTC3_sds = Array(sd_data[:,6:7])

scatter(n_se_tcr[:,1],n_se_tcr[:,2],yerror=n_se_tcr[:,3],markercolor=:black,markersize=5)
scatter!(n_se_1[:,1],n_se_1[:,2],yerror=n_se_1[:,3],markercolor=:blue,markersize=5)
scatter!(n_se_2[:,1],n_se_2[:,2],yerror=n_se_2[:,3],markercolor=:red,markersize=5)
h = scatter!(n_se_3[:,1],n_se_3[:,2],yerror=n_se_3[:,3],markercolor=:cyan,markersize=5,legend=false)

tspan = (0.0,16.0)
t = collect(0.01:0.01:16)

gr()
# bimolecular TCR-pMHC bonds
u₀ = [0.0]
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_T,[])
n_T = Array(solve_model(rates_T))'[2:end]
n_T_sd = n_T .± nT_sds

# TCR-pMHC-CD4 (1)
u₀ = [0.0, 0.0]
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC1,10 .^ rates_T)
n_TC1 = Array(solve_model(rates_TC1))'[2:end,:]
n_TC1_sd = sum(n_TC1 .± nTC1_sds, dims=2)

# TCR-pMHC-CD4 (2)
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC2,10 .^ rates_T)
n_TC2 = Array(solve_model(rates_TC2))'[2:end,:]
n_TC2_sd = sum(n_TC2 .± nTC2_sds, dims=2)

# TCR-pMHC-CD4 (3)
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC3,10 .^ rates_T)
n_TC3 = Array(solve_model(rates_TC3))'[2:end,:]
n_TC3_sd = sum(n_TC3 .± nTC3_sds, dims=2)

plot!(t, n_T_sd, errorbars=:ribbon, alpha=0.1,msc=:black)
plot!(t, n_T, lw=3, lc=:black)
plot!(t,n_TC1_sd, errorbars=:ribbon, alpha=0.1, msc=:blue)
plot!(t, sum(n_TC1,dims=2), lw=3, lc=:blue)
plot!(t,n_TC2_sd, errorbars=:ribbon, alpha=0.1, msc=:red)
plot!(t, sum(n_TC2,dims=2), lw=3, lc=:red)
plot!(t,n_TC3_sd, errorbars=:ribbon, alpha=0.1, msc=:cyan)
plot!(t, sum(n_TC3,dims=2), lw=3, lc=:cyan)

xlabel!("Contact Time tc")
ylabel!("<n>")

#------------------------------------------------------------------------------
# CD4 fitting

cd4_1 = CSV.read("cd4_data_from_muaz.csv")
t_cd41, Pa_cd41, n_se_cd41 = clean_data(cd4_1)
n_cd41 = Pa2n(Pa_cd41)

f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens,cons) # enclose constants

u₀ = [0.0]          # Initial condition
tspan = (0.0,8.0)  # Simulation interval
pinit = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
densities_C = [1700., 850., 1.] #CD4, pMHC, A
cons = []

rates_C,hes_C,l_C,boot_C = optimization(bimolecular!,n_cd41,densities_C',[],model_opt2,t_cd41)

# Setup the ODE problem, then solve
# plotly()
t = collect(0.01:0.01:16)
tspan = (0.0,16.0)
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_C,cons)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_C)
sol_C = Array(solve(prob,Vern7(),saveat=t,abstol=1e-10,reltol=1e-10))'
plot(t,sol_C,lc=:orange,lw=3)
scatter!(n_se_cd41[:,1],n_se_cd41[:,2],yerror=n_se_cd41[:,3],markercolor=:orange)

#-----------------------------------------------------------------------------

densities_T =   [25.0, 38.0, 1.0]       # TCR, pMHC, Ac
densities_TC1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_TC2 = [8.0,  38.0, 17.0, 1.0]
densities_TC3 = [3.0,  38.0, 22.0, 1.0]
densities_C =   [1700., 850., 1.]       # CD4, pMHC, A

n_T_norm = n_T ./25. .*15.
n_TC1_norm = sum(n_TC1,dims=2)
n_TC2_norm = sum(n_TC2,dims=2) ./8. .*15.# ./17. .*10.
n_TC3_norm = sum(n_TC3,dims=2) ./3. .*15.# ./22. .*10.
n_C_norm = sol_C ./1700. .*10. ./850. .*38

plotly()
t = collect(0.01:0.01:16)
plot(t,n_T_norm,lw=3,lc=:black,label="TCR")
plot!(t,n_TC1_norm,lw=3,lc=:blue,label="TCR 15, CD4 10")
plot!(t,n_TC2_norm,lw=3,lc=:red,label="TCR 8, CD4 17")
plot!(t,n_TC3_norm,lw=3,lc=:cyan,label="TCR 3, CD4 22")
plot!(t,n_C_norm,lw=3,lc=:orange,yscale=:linear,label="CD4",legend=:bottomright)
# yaxis!("log<n>")
yaxis!("<n>")
xlabel!("Contact Time tc")
# savefig("n_log_scaled")

coop_indx_TC1 = cooperativity_index.(n_T_norm,n_C_norm,n_TC1_norm)
coop_indx_TC2 = cooperativity_index.(n_T_norm,n_C_norm,n_TC2_norm)
coop_indx_TC3 = cooperativity_index.(n_T_norm,n_C_norm,n_TC3_norm)
plot(t,coop_indx_TC1 .* 100,lw=3,lc=:blue,label="Coop TCR 15, CD4 10")
plot!(t,coop_indx_TC2 .* 100,lw=3,lc=:red,label="Coop TCR 8, CD4 17")
plot!(t,coop_indx_TC3 .* 100,lw=3,lc=:cyan,label="Coop TCR 3, CD4 22")
xlabel!("Contact Time tc")
ylabel!("Cooperative Index (%)")


text:
"""
# Once.jl
"""
module Once

export @once

"""
    @once f()::Type = expression

Define funciton `f()` that evaluates `expression` once on the first call.
Subsequent calls return a cached value without re-evaluating expression.

```
julia> @once f()::String = (println("foo"); "bar")
f (generic function with 1 method)

julia> f()
foo
"bar"

julia> f()
"bar"
```

"""
macro once(ex)
    if !(ex.head == :(=)
    &&   ex.args[1] isa Expr
    &&   ex.args[1].head == :(::)
    &&   ex.args[1].args[1] isa Expr
    &&   ex.args[1].args[1].head == :call
    &&   length(ex.args[1].args[1].args) == 1)
        throw(ArgumentError("Invalid `@once` expression."))
    end

    name = ex.args[1].args[1].args[1]
    type = ex.args[1].args[2]
    init_expression = ex.args[2]

    quote
        let flag = Ref(false), cache = Ref{$type}()

            @noinline init()::$type = $init_expression

            global $name
            Base.@__doc__ @inline function $name()
                if !flag[]
                    cache[] = init()
                    flag[] = true
                end
                cache[]
            end
        end
    end |> esc
end

# Documentation.

readme() = join([
    Docs.doc(@__MODULE__),
    Docs.@doc @once
   ], "\n\n")

end # module


text:
@testset "FairEvaluate Benchmarking" begin
	model = ConstantClassifier()
	X, y = @load_toydata
	results = fairevaluate([model, model], X, y, grp=:Sex, priv_grps=["M"])
	@test length(results)==5
	@test size(results["measures"])==(1,)
	@test size(results["pvalues"])==(1, 2, 2)
	@test size(results["classifier_names"])==(2,)
	@test size(results["results"])==(1, 2, 6)
	@test size(results["tstats"])==(1, 2, 2)
end

function models_fn(fp::FairnessProblem, alpha=1.0)
	########## You can change code segment below to add more models

	rf = @pipeline ContinuousEncoder ConstantClassifier

	models = MLJBase.Model[rf]
	#  Just simply add the base classifiers above and to this list

	model_names = ["RF"] # Add names of models. Preferably keep names short

	################################################################
	model_count = length(models)
	protected_attr = fp.task.grp
	debiasmeasures = fp.task.debiasmeasures
	for i in 1:model_count
    rw = ReweighingSamplingWrapper(classifier=models[i], grp=protected_attr, alpha=alpha)
		lp = LinProgWrapper(classifier=models[i], grp=protected_attr,
								measures=debiasmeasures, alpha=alpha)
		eo = EqOddsWrapper(classifier=models[i], grp=protected_attr, alpha=alpha)
        ce = CalEqOddsWrapper(classifier=models[i], grp=protected_attr, fp_rate=0, fn_rate=1, alpha=alpha)
		append!(models, [rw, lp, eo, ce])
		append!(model_names, model_names[i]*"-".*["Reweighing",
		"LinProg-".*join(debiasmeasures, "-"), "Equalized Odds","Calibrated Equalized Odds"])
	end
	return models, model_names
end

@testset "Get df" begin
	@test size(get_df(student(), models_fn))==(60, 7)
end

@testset "Get Pareto df" begin
	@test size(get_pareto_df(student(), models_fn, 0:1:0.5))==(60, 8)
end


text:
module QueryBuilder

import SearchLight
import DataFrames

import Base: (+) #, select
# TODO: Base.isequal and Base.hash

export from, select, where, limit, offset, order, group, having, prepare

struct MissingModel <: SearchLight.AbstractModel
end

struct QueryPart{T<:SearchLight.AbstractModel}
  model::Type{T}
  query::SearchLight.SQLQuery
end

"""
"""
function from(model::Type{T})::QueryPart{T} where {T<:SearchLight.AbstractModel}
  QueryPart(model, SearchLight.SQLQuery())
end

"""
"""
function select(columns::Vararg{Union{Symbol,String,SearchLight.SQLColumn,SearchLight.SQLRaw}}) :: QueryPart
  QueryPart(MissingModel, SearchLight.SQLQuery(columns = SearchLight.SQLColumns([columns...])))
end

"""
"""
function where(sql_expression::String)::QueryPart
  QueryPart(MissingModel, SearchLight.SQLQuery(where = [SearchLight.SQLWhereExpression(sql_expression)]))
end
function where(sql_expression::String, values::Vararg{Any})::QueryPart
  QueryPart(MissingModel, SearchLight.SQLQuery(where = [SearchLight.SQLWhereExpression(sql_expression, [values...])]))
end

"""
"""
function limit(lim::Int)
  QueryPart(MissingModel, SearchLight.SQLQuery(limit = SearchLight.SQLLimit(lim)))
end

"""
"""
function offset(off::Int)
  QueryPart(MissingModel, SearchLight.SQLQuery(offset = off))
end

"""
"""
function order(ordering::Union{Symbol,String})
  QueryPart(MissingModel, SearchLight.SQLQuery(order = SearchLight.SQLOrder(ordering)))
end
function order(column::Union{Symbol,String}, direction::Union{Symbol,String})
  QueryPart(MissingModel, SearchLight.SQLQuery(order = SearchLight.SQLOrder(column, direction)))
end

"""
"""
function group(columns::Vararg{Union{Symbol,String}})
  QueryPart(MissingModel, SearchLight.SQLQuery(group = SearchLight.SQLColumns([columns...])))
end

"""
"""
function having(sql_expression::String)::QueryPart
  QueryPart(MissingModel, SearchLight.SQLQuery(having = [SearchLight.SQLWhereExpression(sql_expression)]))
end
function having(sql_expression::String, values::Vararg{Any})::QueryPart
  QueryPart(MissingModel, SearchLight.SQLQuery(having = [SearchLight.SQLWhereExpression(sql_expression, [values...])]))
end

"""
"""
function prepare(qb::QueryPart{T}) where {T<:SearchLight.AbstractModel}
  (qb.model::Type{T}, qb)
end
function prepare(model::Type{T}, qb::QueryPart) where {T<:SearchLight.AbstractModel}
  prepare(from(model) + qb)
end

"""
"""
function (+)(q::SearchLight.SQLQuery, r::SearchLight.SQLQuery)
  SearchLight.SQLQuery(
    columns = vcat(q.columns, r.columns),
    where   = vcat(q.where, r.where),
    limit   = r.limit.value == SearchLight.SQLLimit_ALL ? q.limit : r.limit,
    offset  = r.offset != 0 ? r.offset : q.offset,
    order   = vcat(q.order, r.order),
    group   = vcat(q.group, r.group),
    having  = vcat(q.having, r.having)
  )
end

"""
"""
function (+)(q::QueryPart, r::QueryPart)
  QueryPart(
    r.model == MissingModel ? q.model : r.model,
    q.query + r.query
  )
end

### API

"""
"""
function DataFrames.DataFrame(m::Type{T}, qp::QueryBuilder.QueryPart, j::Union{Nothing,Vector{SearchLight.SQLJoin{N}}} = nothing)::DataFrames.DataFrame where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  SearchLight.DataFrame(m, qp.query, j)
end

"""
"""
function SearchLight.find(m::Type{T}, qp::QueryBuilder.QueryPart,
                      j::Union{Nothing,Vector{SearchLight.SQLJoin{N}}} = nothing)::Vector{T} where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  SearchLight.find(m, qp.query, j)
end

function SearchLight.first(m::Type{T}, qp::QueryBuilder.QueryPart)::Union{Nothing,T} where {T<:SearchLight.AbstractModel}
  SearchLight.find(m, qp + QueryBuilder.limit(1)) |> onereduce
end

function SearchLight.last(m::Type{T}, qp::QueryBuilder.QueryPart)::Union{Nothing,T} where {T<:SearchLight.AbstractModel}
  SearchLight.find(m, qp + QueryBuilder.limit(1)) |> onereduce
end

function SearchLight.count(m::Type{T}, qp::QueryBuilder.QueryPart)::Int where {T<:SearchLight.AbstractModel}
  SearchLight.count(m, qp.query)
end

function SearchLight.sql(m::Type{T}, qp::QueryBuilder.QueryPart)::String where {T<:SearchLight.AbstractModel}
  SearchLight.sql(m, qp.query)
end

end


text:
function mask!(a,pad)
    for i in 1:size(a,1)
        for j in size(a,2):-1:1
            if a[i,j] != pad
                break
            elseif j > 1 && a[i,j-1] == pad
                a[i,j] = 0
            end
        end
    end
    return a
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
#=

@author : Spencer Lyon <spencer.lyon@nyu.edu>

=#
using QuantEcon
using Plots
pyplot()

# model parameters
a_0 = 100
a_1 = 0.5
ρ = 0.9
σ_d = 0.05
β = 0.95
c = 2
γ = 50.0
θ = 0.002
ac = (a_0 - c) / 2.0

# Define LQ matrices
R = [ 0   ac    0;
     ac -a_1  0.5;
     0.  0.5    0]
R = -R  # For minimization
Q = Matrix([γ / 2.0]')
A = [1. 0. 0.;
     0. 1. 0.;
     0. 0.  ρ]
B = [0. 1. 0.]'
C = [0. 0. σ_d]'

## Functions

function evaluate_policy(θ::AbstractFloat, F::AbstractArray)
    rlq = RBLQ(Q, R, A, B, C, β, θ)
    K_F, P_F, d_F, O_F, o_F = evaluate_F(rlq, F)
    x0 = [1.0 0.0 0.0]'
    value = - x0' * P_F * x0 - d_F
    entropy = x0' * O_F * x0 + o_F
    return value[1], entropy[1]  # return scalars
end

function value_and_entropy{TF<:AbstractFloat}(emax::AbstractFloat,
                                              F::AbstractArray{TF},
                                              bw::String,
                                              grid_size::Integer=1000)
    if lowercase(bw) == "worst"
        θs = 1 ./ linspace(1e-8, 1000, grid_size)
    else
        θs = -1 ./ linspace(1e-8, 1000, grid_size)
    end

    data = Array{TF}(grid_size, 2)

    for (i, θ) in enumerate(θs)
        data[i, :] = collect(evaluate_policy(θ, F))
        if data[i, 2] >= emax  # stop at this entropy level
            data = data[1:i, :]
            break
        end
    end
    return data
end

## Main

# compute optimal rule
optimal_lq = LQ(Q, R, A, B, C, zero(B'A), bet=β)
Po, Fo, Do = stationary_values(optimal_lq)

# compute robust rule for our θ
baseline_robust = RBLQ(Q, R, A, B, C, β, θ)
Fb, Kb, Pb = robust_rule(baseline_robust)

# Check the positive definiteness of worst-case covariance matrix to
# ensure that θ exceeds the breakdown point
test_matrix = eye(size(Pb, 1)) - (C' * Pb * C ./ θ)[1]
eigenvals, eigenvecs = eig(test_matrix)
@assert all(eigenvals .>= 0)

emax = 1.6e6

# compute values and entropies
optimal_best_case = value_and_entropy(emax, Fo, "best")
robust_best_case = value_and_entropy(emax, Fb, "best")
optimal_worst_case = value_and_entropy(emax, Fo, "worst")
robust_worst_case = value_and_entropy(emax, Fb, "worst")

# we reverse order of "worst_case"s so values are ascending
data_pairs = ((optimal_best_case, optimal_worst_case),
              (robust_best_case, robust_worst_case))

egrid = linspace(0, emax, 100)
egrid_data = Array{Float64}[]
for data_pair in data_pairs
    for data in data_pair
        x, y = data[:, 2], data[:, 1]
        curve = LinInterp(x, y)
        push!(egrid_data, curve.(egrid))
    end
end
plot(egrid, egrid_data, color=[:red :red :blue :blue])
plot!(egrid, egrid_data[1], fillrange=egrid_data[2],
      fillcolor=:red, fillalpha=0.1, color=:red, legend=:none)
plot!(egrid, egrid_data[3], fillrange=egrid_data[4],
      fillcolor=:blue, fillalpha=0.1, color=:blue, legend=:none)
plot!(xlabel="Entropy", ylabel="Value")


text:
# This file was generated, do not modify it. # hide
boston.Crim_x_Zn = boston.Crim .* boston.Zn;