text:
using Test
using Random: seed!
using StaticArrays: SVector
using QuantumLattices.Essentials.Spatials
using QuantumLattices.Essentials: dtype, kind, reset!
using QuantumLattices.Essentials.QuantumNumbers: Momentum₁, Momentum₂, Momentum₃, AbelianNumbers
using QuantumLattices.Interfaces: decompose, rank, dimension, expand
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Prerequisites.Traits: contentnames, getcontent
using QuantumLattices.Prerequisites.SimpleTrees: leaves

@testset "distance" begin
    @test distance([0.0, 0.0], [1.0, 1.0]) ≈ sqrt(2.0)
end

@testset "azimuthd" begin
    @test azimuthd([+1.0]) ≈ 0.0
    @test azimuthd([-1.0]) ≈ 180.0
    @test azimuthd([+1.0, +1.0]) ≈ 45.0
    @test azimuthd([-1.0, -1.0]) ≈ 225.0
    @test azimuthd([+1.0, +1.0, -2.0]) ≈ 45.0
    @test azimuthd([+1.0, -1.0, +2.0]) ≈ 315.0
end

@testset "azimuth" begin
    @test azimuth([+1.0]) ≈ 0.0
    @test azimuth([-1.0]) ≈ pi
    @test azimuth([+1.0, +1.0]) ≈ 1//4*pi
    @test azimuth([-1.0, -1.0]) ≈ 5//4*pi
    @test azimuth([+1.0, +1.0, -2.0]) ≈ 1//4*pi
    @test azimuth([+1.0, -1.0, +2.0]) ≈ 7//4*pi
end

@testset "polard" begin
    @test polard([1.0, 0.0, 1.0]) ≈ 45.0
end

@testset "polard" begin
    @test polar([1.0, 0.0, 1.0]) ≈ 1//4*pi
end

@testset "volume" begin
    @test volume([[1]]) == volume([1]) == volume([1, 0]) == volume([1, 0, 0]) == 1
    @test volume([1], [2]) == 0
    @test volume([[1, 0], [0, 1]]) == volume([1, 0], [0, 1]) == volume([1, 0, 0], [0, 1, 0]) == 1
    @test volume([1], [2], [3]) == volume([1, 0], [0, 1], [1, 1]) == 0
    @test volume([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) == volume([1, 0, 0], [0, 1, 0], [0, 0, 1]) == 1
end

@testset "isparallel" begin
    @test isparallel([0.0, 0.0], [1.0, 1.0]) == 1
    @test isparallel([2.0, 2.0], [1.0, 1.0]) == 1
    @test isparallel([-2.0, -2.0], [1.0, 1.0]) == -1
    @test isparallel([1.0, -1.0], [1.0, 1.0]) == 0
end

@testset "isonline" begin
    seed!()
    p1, p2 = [0.0, 0.0], [1.0, 1.0]
    @test isonline([0.0, 0.0], p1, p2, ends = (true, rand(Bool))) == true
    @test isonline([0.0, 0.0], p1, p2, ends = (false, rand(Bool))) == false
    @test isonline([1.0, 1.0], p1, p2, ends = (rand(Bool), true)) == true
    @test isonline([1.0, 1.0], p1, p2, ends = (rand(Bool), false)) == false
    @test isonline([0.5, 0.5], p1, p2, ends = (rand(Bool), rand(Bool))) == true
    @test isonline([1.1, 1.1], p1, p2, ends = (rand(Bool), rand(Bool))) == false
end

@testset "decompose" begin
    a, c = randn(3), randn()
    @test all(decompose(c*a, a) .≈ (c,))

    a1, a2 = randn(2), randn(2)
    c1, c2 = randn(2)
    @test all(decompose(c1*a1+c2*a2, a1, a2) .≈ (c1, c2))

    a1, a2 = randn(3), randn(3)
    c1, c2 = randn(2)
    @test all(decompose(c1*a1+c2*a2, a1, a2) .≈ (c1, c2))

    a1, a2, a3 = randn(3), randn(3), randn(3)
    c1, c2, c3 = randn(3)
    @test all(decompose(c1*a1+c2*a2+c3*a3, a1, a2, a3) .≈ (c1, c2, c3))
end

@testset "isintratriangle" begin
    seed!()
    p1, p2, p3 = [-1.0, 0.0], [0.0, 1.0], [1.0, 0.0]
    for (i, p) in enumerate((p1, p2, p3))
        @test isintratriangle(p, p1, p2, p3, vertexes=Tuple(j == i ? true : rand(Bool) for j = 1:3), edges=(rand(Bool), rand(Bool), rand(Bool))) == true
        @test isintratriangle(p, p1, p2, p3, vertexes=Tuple(j == i ? false : rand(Bool) for j = 1:3), edges=(rand(Bool), rand(Bool), rand(Bool))) == false
    end
    for (i, p) in enumerate(((p1+p2)/2, (p2+p3)/2, (p3+p1)/2))
        @test isintratriangle(p, p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=Tuple(j == i ? true : rand(Bool) for j = 1:3)) == true
        @test isintratriangle(p, p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=Tuple(j == i ? false : rand(Bool) for j = 1:3)) == false
    end
    @assert isintratriangle([0.0, 0.5], p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=(rand(Bool), rand(Bool), rand(Bool))) == true
    @assert isintratriangle([0.0, 1.5], p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=(rand(Bool), rand(Bool), rand(Bool))) == false
end

@testset "issubordinate" begin
    @test issubordinate([2.0], [[1.0]]) == true
    @test issubordinate([1.5], [[1.0]]) == false
    @test issubordinate([0.5, 3.5], [[0.5, 0.5], [-0.5, 0.5]]) == true
    @test issubordinate([0.5, 3.0], [[0.5, 0.5], [-0.5, 0.5]]) == false
    @test issubordinate([4.0, 4.0, 4.0], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) == true
    @test issubordinate([4.0, 4.0, 2.5], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) == false
end

@testset "reciprocals" begin
    @test reciprocals([[1.0]]) ≈ [[1.0]]*2pi
    @test reciprocals([[1.0, 0.0], [0.0, 1.0]]) ≈ [[1.0, 0.0], [0.0, 1.0]]*2pi
    @test reciprocals([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) ≈ [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]*2pi
end

@testset "translate" begin
    @test translate([1.0 2.0; 3.0 4.0], [1.0, 2.0]) == [2.0 3.0; 5.0 6.0]
end

@testset "rotate" begin
    @test rotate([1.0 2.0; 1.0 2.0], pi/4, axis=([1.0, 1.0], (0.0, 0.0))) ≈ [1.0 1.0; 1.0 1.0+√2.0]
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([0.0, 1.0, 1.0], (pi/2, 0.0))) ≈ reshape([2.0, 1.0+√2/2, 1.0+√6/2], 3, 1)
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([1.0, 0.0, 1.0], (pi/2, pi/2))) ≈ reshape([1.0+√6/2, 2.0, 1.0+√2/2], 3, 1)
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([1.0, 1.0, 0.0], (0.0, 0.0))) ≈ reshape([1.0+√2/2, 1.0+√6/2, 2.0], 3, 1)
end

@testset "tile" begin
    @test tile(reshape([0.0, 0.0], 2, 1), [[1.0, 0.0], [0.0, 1.0]], ((0.5, 0.5), (-0.5, -0.5))) == [0.5 -0.5; 0.5 -0.5]
end

@testset "minimumlengths" begin
    @test minimumlengths(reshape([0.0, 0.0], 2, 1), [[1.0, 0.0], [0.0, 1.0]], 7) ≈ [1.0, √2, 2.0, √5, 2*√2, 3.0, √10]
end

@testset "intralinks" begin
    ps, a1, a2 = [0.0 0.5; 0.0 0.5], [1.0, 0.0], [0.0, 1.0]
    neighbors = Dict(1=>1.0)
    links = [(1, 1, 1, SVector(0.0, -1.0)), (1, 1, 1, SVector(-1.0, 0.0)), (1, 2, 2, SVector(0.0, -1.0)), (1, 2, 2, SVector(-1.0, 0.0))]
    @test intralinks(ps, [a1, a2], neighbors) == links
end

@testset "interlinks" begin
    ps1 = [0.0 1.0; 0.0 0.0]
    ps2 = [0.0 1.0; 1.0 1.0]
    neighbors = Dict(1=>1.0)
    links = [(1, 1, 1), (1, 2, 2)]
    @test interlinks(ps1, ps2, neighbors) == links
end

@testset "CPID" begin
    @test CPID(1) == CPID('T', 1) == CPID(scope='T', site=1)
    @test CPID(scope="tz", site=1)|>string == "CPID(\"tz\", 1)"
end

@testset "Point" begin
    point = Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))
    @test point|>dtype == point|>typeof|>dtype == Float
    @test point == Point(CPID(1, 1), [0.0, 0.0], [0.0, 0.0])
    @test isequal(point, Point{2}(CPID(1, 1), [0.0, 0.0], [0.0, 0.0]))
    @test point|>length == point|>typeof|>length == 1
    @test point|>eltype == point|>typeof|>eltype == Point{2, CPID{Int}, Float}
    @test point|>rank == point|>typeof|>rank == 1
    @test point|>pidtype == point|>typeof|>pidtype == CPID{Int}
    @test point|>dimension == point|>typeof|>dimension == 2
    @test point|>kind == point|>typeof|>kind == 0
    @test point|>string == "Point(CPID(1, 1), [0.0, 0.0], [0.0, 0.0])"
    @test point[1] == point
    @test point|>collect == [point]

    @test isintracell(Point(PID(1), (0.0, 0.0), (0.0, 0.0))) == true
    @test isintracell(Point(PID(1), (0.0, 0.0), (1.0, 0.0))) == false
end

@testset "Bond" begin
    bond = Bond(1, Point(PID(1), (0.0, 0.0), (0.0, 0.0)), Point(PID(2), (0.0, 1.0), (0.0, 1.0)))
    @test bond|>deepcopy == bond
    @test isequal(bond|>deepcopy, bond)
    @test bond|>string == "Bond(1, Point(PID(1), [0.0, 0.0], [0.0, 0.0]), Point(PID(2), [0.0, 1.0], [0.0, 1.0]))"
    @test bond|>reverse == Bond(1, Point(PID(2), (0.0, 1.0), (0.0, 1.0)), Point(PID(1), (0.0, 0.0), (0.0, 0.0)))
    @test bond|>length == bond|>typeof|>length == 2
    @test bond|>eltype == bond|>typeof|>eltype == Point{2, PID, Float}
    @test bond|>rank == bond|>typeof|>rank == 2
    @test bond|>pidtype == bond|>typeof|>pidtype == PID
    @test bond|>dimension == bond|>typeof|>dimension == 2
    @test bond|>kind == 1
    @test bond|>rcoord == [0.0, 1.0]
    @test bond|>icoord == [0.0, 1.0]
    @test bond|>isintracell == false
    @test bond[1] == bond.epoint && bond[2] == bond.spoint
    @test bond|>collect == [Point(PID(2), (0.0, 1.0), (0.0, 1.0)), Point(PID(1), (0.0, 0.0), (0.0, 0.0))]
end

@testset "Lattice" begin
    lattice = Lattice("Tuanzi", [Point(PID(1), (0.5, 0.5), (0.0, 0.0))], vectors=[[1.0, 0.0], [0.0, 1.0]], neighbors=1)
    @test lattice|>typeof|>contentnames == (:name, :pids, :rcoords, :icoords, :vectors, :reciprocals, :neighbors)
    @test lattice|>deepcopy == lattice
    @test isequal(lattice|>deepcopy, lattice)
    @test lattice|>string == "Lattice(Tuanzi)\n  with 1 point:\n    Point(PID(1), [0.5, 0.5], [0.0, 0.0])\n  with 2 translation vectors:\n    [1.0, 0.0]\n    [0.0, 1.0]\n  with 1 order of nearest neighbors:\n    1 => 1.0\n"
    @test lattice|>length == 1
    @test lattice|>dimension == lattice|>typeof|>dimension == 2
    @test lattice|>keytype == lattice|>typeof|>keytype == PID
    @test lattice|>valtype == lattice|>typeof|>valtype == Point{2, PID, Float}
    @test lattice|>dtype == lattice|>typeof|>dtype == Float
    @test lattice[LatticeIndex{'R'}(1)] == lattice[LatticeIndex{'R'}(PID(1))] == SVector(0.5, 0.5)
    @test lattice[LatticeIndex{'I'}(1)] == lattice[LatticeIndex{'I'}(PID(1))] == SVector(0.0, 0.0)
    @test lattice[LatticeIndex{'P'}(1)] == lattice[LatticeIndex{'P'}(PID(1))] == Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    @test lattice|>nneighbor == 1

    zerothbs = [Point(PID(1), (0.5, 0.5), (0.0, 0.0))]
    insidebs = []
    acrossbs = [Bond(1, Point(PID(1), [0.5, 0.5], [0.0, 0.0]), Point(PID(1), [0.5, -0.5], [0.0, -1.0])),
                Bond(1, Point(PID(1), [0.5, 0.5], [0.0, 0.0]), Point(PID(1), [-0.5, 0.5], [-1.0, 0.0]))
    ]
    @test bonds(lattice, zerothbonds) == zerothbs
    @test bonds(lattice, insidebonds) == insidebs
    @test bonds(lattice, acrossbonds) == acrossbs
    @test setdiff(bonds(lattice), [zerothbs; insidebs; acrossbs])|>length == 0

    lattice = Lattice("SuperTuanzi",
                [Lattice("Tuanzi1", [Point(CPID(1, 1), (0.0, 0.0))], neighbors=0),
                Lattice("Tuanzi2", [Point(CPID(2, 1), (0.5, 0.5))], neighbors=0),
                ],
                vectors=[[1.0, 0.0], [0.0, 1.0]],
                neighbors=2,
                )
    @test setdiff(bonds(lattice), bonds(lattice, zerothbonds, insidebonds, acrossbonds))|>length == 0
end

@testset "SuperLattice" begin
    lattice = SuperLattice("SuperTuanzi",
                [Lattice("TuanziSys", [Point(CPID(1, 1), (0.0, 0.0)), Point(CPID(1, 2), (0.5, 0.5))], neighbors=Dict(1=>√2/2)),
                Lattice("TuanziEnv", [Point(CPID(2, 1), (-0.05, -0.05)), Point(CPID(2, 2), (0.55, 0.55))], neighbors=Dict{Int, Float}()),
                ],
                vectors=[[1.0, 0.0], [0.0, 1.0]],
                neighbors=Dict(1=>√2/2, -1=>√2/20)
                )
    @test lattice|>typeof|>contentnames == (:sublattices, :name, :pids, :rcoords, :icoords, :vectors, :reciprocals, :neighbors)
    @test latticetype(lattice) == latticetype(typeof(lattice)) == Lattice{2, CPID{Int}, Float}
    @test setdiff(bonds(lattice), bonds(lattice, zerothbonds, insidebonds, acrossbonds))|>length == 0
    @test setdiff(bonds(lattice, insidebonds), bonds(lattice, intrabonds, interbonds))|>length == 0
end

@testset "Cylinder" begin
    cylinder = Cylinder{CPID{String}}("Tuanzi", [0.0 0.0; 0.0 1.0], SVector(1.0, 0.0), vector=[0.0, 2.0], neighbors=1)
    @test cylinder|>typeof|>contentnames == (:block, :translation, :name, :pids, :rcoords, :icoords, :vectors, :reciprocals, :neighbors)
    insert!(cylinder, "A", "B")
    insert!(cylinder, "C3", cut=2, scopes=["C1", "C1", "C2", "C2"])
    @test cylinder.pids == [CPID("C1", 1), CPID("C3", 2), CPID("C3", 1), CPID("C1", 2), CPID("C2", 1), CPID("C2", 2)]
    @test cylinder.rcoords == [-1.0 -1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0]
    @test cylinder.icoords == [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]
    lattice = cylinder("C1", "C2", "C3")
    @test lattice.pids == [CPID("C1", 1), CPID("C1", 2), CPID("C2", 1), CPID("C2", 2), CPID("C3", 1), CPID("C3", 2)]
    @test lattice.rcoords == [-1.0 -1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0]
    @test lattice.icoords == [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]
end

@testset "LatticeBonds" begin
    @test eltype(typeof(allbonds)) === AbstractBond
    @test eltype(typeof(zerothbonds)) === Point
    @test eltype(typeof(insidebonds)) === eltype(typeof(acrossbonds)) === Bond
    @test eltype(typeof(intrabonds)) === eltype(typeof(interbonds)) === Bond
    @test eltype(AbstractLattice{2, PID, Float}, allbonds) == eltype(AbstractLattice{2, PID, Float}, Val(allbonds)) == AbstractBond{2, PID, Float}
    @test eltype(Lattice{3, CPID{Char}, Float}, zerothbonds) == eltype(Lattice{3, CPID{Char}, Float}, Val(zerothbonds)) == Point{3, CPID{Char}, Float}
    @test eltype(Lattice{1, CPID{String}, Float}, insidebonds) == eltype(Lattice{1, CPID{String}, Float}, Val(insidebonds)) == Bond{1, CPID{String}, Float}

    @test expand(Lattice, Val((allbonds,))) == Lattice|>latticebondsstructure|>leaves|>Tuple
    @test expand(Lattice, Val((zerothbonds,))) == (zerothbonds,)
    @test expand(Lattice, Val((insidebonds,))) == (insidebonds,)
    @test expand(Lattice, Val((acrossbonds,))) == (acrossbonds,)

    @test expand(SuperLattice, Val((allbonds,))) == SuperLattice|>latticebondsstructure|>leaves|>Tuple
    @test expand(SuperLattice, Val((zerothbonds,))) == (zerothbonds,)
    @test expand(SuperLattice, Val((insidebonds,))) == (intrabonds, interbonds)
    @test expand(SuperLattice, Val((acrossbonds,))) == (acrossbonds,)
    @test expand(SuperLattice, Val((intrabonds,))) == (intrabonds,)
    @test expand(SuperLattice, Val((interbonds,))) == (interbonds,)
end

@testset "Bonds" begin
    lattice = Lattice("Tuanzi", [Point(PID(1), (0.5, 0.5), (0.0, 0.0))], vectors=[[1.0, 0.0], [0.0, 1.0]], neighbors=1)
    bs = Bonds(lattice)
    @test bs|>eltype == bs|>typeof|>eltype == AbstractBond{2, PID, Float}
    @test bs==Bonds(lattice, allbonds) && isequal(bs, Bonds(lattice, allbonds))
    @test (bs|>length == 3) && (bs|>size == (3,))
    @test summary(bs) == "3-element Bonds"
    @test bs|>bondtypes == bs|>typeof|>bondtypes == (zerothbonds, insidebonds, acrossbonds)
    @test bs|>latticetype == bs|>typeof|>latticetype == Lattice{2, PID, Float}
    @test bs|>rank == bs|>typeof|>rank == 3
    @test bs|>collect == bonds(lattice)
    @test (bs[1] == bs.bonds[1][1]) && (bs[2] == bs.bonds[3][1]) && (bs[3] == bs.bonds[3][2])
    @test filter(allbonds, bs, :include|>Val) == filter((zerothbonds, insidebonds, acrossbonds), bs, :include|>Val) == bs
    @test filter(zerothbonds, bs, :include|>Val) == filter((insidebonds, acrossbonds), bs, :exclude|>Val) == Bonds(lattice, zerothbonds)
    @test filter(insidebonds, bs, :include|>Val) == filter((zerothbonds, acrossbonds), bs, :exclude|>Val) == Bonds(lattice, insidebonds)
    @test filter(acrossbonds, bs, :include|>Val) == filter((zerothbonds, insidebonds), bs, :exclude|>Val) == Bonds(lattice, acrossbonds)
    emptybs = Bonds{(zerothbonds, insidebonds, acrossbonds), Lattice{2, PID}}((Point{2, PID, Float}[], Bond{2, PID, Float}[], Bond{2, PID, Float}[]))
    @test filter(bond->(dimension(bond) == 3), bs) == empty!(deepcopy(bs)) == empty(bs) == emptybs
    @test reset!(emptybs, lattice) == bs
end

@testset "Segment" begin
    segment = Segment(SVector(0.0, 0.0), SVector(1.0, 1.0), 10)
    @test segment == Segment((0.0, 0.0), (1.0, 1.0), 10)
    @test isequal(segment,  Segment((0.0, 0.0), (1.0, 1.0), 10))
    @test size(segment) == (10,)
    @test string(segment) == "[p₁, p₂) with p₁=[0.0, 0.0] and p₂=[1.0, 1.0]"

    segment = Segment(1, 5, 5, ends=(true, true))
    @test string(segment) == "[1.0, 5.0]"
    @test segment[1]==1 && segment[end]==5
    for (i, seg) in enumerate(segment)
        @test seg ≈ segment[i]
    end

    segment = Segment(1, 5, 4, ends=(true, false))
    @test string(segment) == "[1.0, 5.0)"
    @test segment[1]==1 && segment[end]==4
    for (i, seg) in enumerate(segment)
        @test seg ≈ segment[i]
    end

    segment = Segment(1, 5, 4, ends=(false, true))
    @test string(segment) == "(1.0, 5.0]"
    @test segment[1]==2 && segment[end]==5
    for (i, seg) in enumerate(segment)
        @test seg ≈ segment[i]
    end

    segment = Segment(1, 5, 3, ends=(false, false))
    @test string(segment) == "(1.0, 5.0)"
    @test segment[1]==2 && segment[end]==4
    for (i, seg) in enumerate(segment)
        @test seg ≈ segment[i]
    end
end

@testset "BrillouinZone" begin
    @test contentnames(BrillouinZone) == (:reciprocals, :contents)

    recipls = reciprocals([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
    momenta = AbelianNumbers('C', [Momentum₂{9, 10}(j-1, i-1) for (i, j) in Iterators.flatten((Iterators.product(1:10, 1:9),))], collect(0:90), :indptr)
    bz = BrillouinZone(recipls, momenta)
    @test getcontent(bz, Val(:contents)) == (bz.momenta,)
    @test dtype(bz) == dtype(typeof(bz)) == Float
    @test bz == BrillouinZone(Momentum₂{9, 10}, recipls)

    recipls = reciprocals([[1.0, 0.0, 0.0]])
    @test BrillouinZone{:q}(Momentum₁{10}, recipls) == BrillouinZone{:q}(recipls, 10)

    recipls = reciprocals([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
    @test BrillouinZone(Momentum₂{10, 10}, recipls) == BrillouinZone(recipls, 10)

    recipls = reciprocals([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
    @test BrillouinZone(Momentum₃{10, 10, 10}, recipls) == BrillouinZone(recipls, 10)
end

@testset "ReciprocalZone" begin
    @test contentnames(ReciprocalZone) == (:contents, :volume)

    rz = ReciprocalZone([[1.0]], len=10)
    @test getcontent(rz, :contents) == (rz.momenta,)
    @test rz == ReciprocalZone(rz.momenta, rz.volume)
    @test rz == ReciprocalZone([[1.0]], Segment(-1//2, 1//2, 10))
    @test rz == ReciprocalZone([[1.0]], (Segment(-1//2, 1//2, 10),))

    rz = ReciprocalZone{:q}([[1.0]], len=10)
    @test rz == ReciprocalZone{:q}(rz.momenta, rz.volume)
    @test rz == ReciprocalZone{:q}([[1.0]], Segment(-1//2, 1//2, 10))
    @test rz == ReciprocalZone{:q}([[1.0]], (Segment(-1//2, 1//2, 10),))

    b₁, b₂, b₃ = [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]
    @test ReciprocalZone([b₁], Segment(-1, +1, 10)).volume == 2
    @test ReciprocalZone([b₁, b₂], Segment(-1, +1, 10), Segment(-1, +1, 10)).volume == 4
    @test ReciprocalZone([b₁, b₂, b₃], Segment(-1, +1, 10), Segment(-1, +1, 10), Segment(-1, +1, 10)).volume == 8

    rz = ReciprocalZone(BrillouinZone{:q}(Momentum₂{8, 8}, [b₁, b₂]))
    @test rz == ReciprocalZone{:q}([b₁, b₂], Segment(0, 1, 8), Segment(0, 1, 8))
end

@testset "ReciprocalPath" begin
    @test contentnames(ReciprocalPath) == (:contents,)

    b₁, b₂ = [1.0, 0.0], [0.0, 1.0]
    s₁ = Segment((0.0, 0.0), (0.5, 0.0), 100)
    s₂ = Segment((0.5, 0.0), (0.5, 0.5), 100)
    s₃ = Segment((0.5, 0.5), (0.0, 0.0), 100)

    rp = ReciprocalPath([b₁, b₂], s₁, s₂, s₃)
    @test getcontent(rp, :contents) == (rp.momenta,)
    @test rp == ReciprocalPath(rp.momenta)
    @test rp == ReciprocalPath([b₁, b₂], (s₁, s₂, s₃))

    rp = ReciprocalPath{:q}([b₁, b₂], s₁, s₂, s₃)
    @test rp == ReciprocalPath{:q}(rp.momenta)
    @test rp == ReciprocalPath{:q}([b₁, b₂], (s₁, s₂, s₃))

    rp = ReciprocalPath([b₁+b₂], line"X₂-X₁", len=10)
    @test rp == ReciprocalPath([b₁+b₂], (-1//2,)=>(1//2,), len=10)

    rp = ReciprocalPath([b₁, b₂], rectangle"Γ-X-M-Γ", len=10)
    @test rp == ReciprocalPath([b₁, b₂], (0//1, 0//1)=>(1//2, 0//1), (1//2, 0//1)=>(1//2, 1//2), (1//2, 1//2)=>(0//1, 0//1), len=10)

    rp = ReciprocalPath{:q}([b₁, b₂], hexagon"Γ-K-M-Γ", len=10)
    @test rp == ReciprocalPath{:q}([b₁, b₂], (0//1, 0//1)=>(1//3, 1//3), (1//3, 1//3)=>(1//2, 0//1), (1//2, 0//1)=>(0//1, 0//1), len=10)
end


text:
"""
Interface to libpq - which interfaces with PostgreSQL's backend server

All functions should be considered unsafe (will segfault with bad pointers.)
Also, pointers's need their memory freed by calling the right PQ* functions.
"""
macro c(ret_type, func, arg_types, lib)
    local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
    quote
        $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ),
                                            $ret_type, $arg_types, $(args_in...) )
    end
end

abstract PGconn
abstract PGresult
abstract PGcancel

typealias ResultPtr Nullable{Ptr{Libpq.PGresult}}
typealias ConnPtr Nullable{Ptr{Libpq.PGconn}}

typealias Oid UInt32
typealias ConnStatusType UInt32

const connection_status = Dict(
    (const CONNECTION_OK = 0)                   => :ok,
    (const CONNECTION_BAD = 1)                  => :bad,
    (const CONNECTION_STARTED = 2)              => :started,
    (const CONNECTION_MADE = 3)                 => :made,
    (const CONNECTION_AWAITING_RESPONSE = 4)    => :awaiting_response,
    (const CONNECTION_AUTH_OK = 5)              => :auth_ok,
    (const CONNECTION_SETENV = 6)               => :setenv,
    (const CONNECTION_SSL_STARTUP = 7)          => :ssl_startup,
    (const CONNECTION_NEEDED = 8)               => :needed,
)
typealias ExecStatusType UInt32
const exec_status = Dict(
    (const PGRES_EMPTY_QUERY = 0)               => :empty_query,
    (const PGRES_COMMAND_OK = 1)                => :command_ok,
    (const PGRES_TUPLES_OK = 2)                 => :tuples_ok,
    (const PGRES_COPY_OUT = 3)                  => :copy_out,
    (const PGRES_COPY_IN = 4)                   => :copy_in,
    (const PGRES_BAD_RESPONSE = 5)              => :bad_response,
    (const PGRES_NONFATAL_ERROR = 6)            => :nonfatal_error,
    (const PGRES_FATAL_ERROR = 7)               => :fatal_error,
    (const PGRES_COPY_BOTH = 8)                 => :copy_both,
    (const PGRES_SINGLE_TUPLE = 9)              => :single_tuple,
)

const error_field = Dict(
    :severity		            => 'S', # ERROR, WARNING etc.
    :sqlstate		            => 'C', # see apendix A in the pg man
    :message_primary            => 'M',
    :message_detail	            => 'D',
    :message_hint	            => 'H',
    :statement_position         => 'P',
    :internal_position          => 'p',
    :internal_query	            => 'q',
    :context			        => 'W',
    :schema_name		        => 's',
    :table_name		            => 't',
    :column_name		        => 'c',
    :datatype_name	            => 'd',
    :constraint_name            => 'n',
    :source_file		        => 'F',
    :source_line		        => 'L',
    :source_function            => 'R',
)

const error_state = Dict(
    "00" => :successful_completion,
    "01" => :warning,
    "02" => :no_data,
    "03" => :sql_statement_not_yet_complete,
    "08" => :connection_exception,
    "09" => :triggered_action_exception,
    "0A" => :feature_not_supported,
    "0B" => :invalid_transaction_initiation,
    "0F" => :locator_exception,
    "0L" => :invalid_grantor,
    "0P" => :invalid_role_specification,
    "0Z" => :diagnostics_exception,
    "20" => :case_not_found,
    "21" => :cardinality_violation,
    "22" => :data_exception,
    "23" => :integrity_constraint_violation,
    "24" => :invalid_cursor_state,
    "25" => :invalid_transaction_state,
    "26" => :invalid_sql_statement_name,
    "27" => :triggered_data_change_violation,
    "28" => :invalid_authorization_specification,
    "2B" => :dependent_privilege_descriptors_still_exist,
    "2D" => :invalid_transaction_termination,
    "2F" => :sql_routine_exception,
    "34" => :invalid_cursor_name,
    "38" => :external_routine_exception,
    "39" => :external_routine_invocation_exception,
    "3B" => :savepoint_exception,
    "3D" => :invalid_catalog_name,
    "3F" => :invalid_schema_name,
    "40" => :transaction_rollback,
    "42" => :syntax_error_or_access_rule_violation,
    "44" => :with_check_option_violation,
    "53" => :insufficient_resources,
    "54" => :program_limit_exceeded,
    "55" => :object_not_in_prerequisite_state,
    "57" => :operator_intervention,
    "58" => :system_error,
    "F0" => :configuration_file_error,
    "HV" => :foreign_data_wrapper_error,
    "P0" => :plpgsql_error,
    "XX" => :internal_error,
)

#### CONNECTIONS
@c Ptr{PGconn}      PQsetdbLogin (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8},
                              Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}) libpq
@c Ptr{PGconn}      PQconnectdb (Ptr{UInt8},) libpq
@c Void             PQfinish (Ptr{PGconn},) libpq
@c Ptr{UInt8}       PQerrorMessage (Ptr{PGconn},) libpq
@c ExecStatusType   PQresultStatus (Ptr{PGresult},) libpq
@c ConnStatusType   PQstatus (Ptr{PGconn},) libpq
@c Ptr{UInt8}       PQresultErrorMessage (Ptr{PGresult},) libpq
@c Ptr{UInt8}       PQresultErrorField (Ptr{PGresult}, Cint) libpq
@c Void             PQsetNoticeReceiver (Ptr{PGconn}, Ptr{Void}, Ptr{Void}) libpq

#### EXEC COMMANDS
@c Ptr{PGresult}    PQexec (Ptr{PGconn}, Ptr{UInt8}) libpq
@c Ptr{PGresult}    PQgetResult (Ptr{PGconn},) libpq # for the end of a copy command
@c Cint             PQputCopyData (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
@c Cint             PQputCopyEnd (Ptr{PGconn}, Ptr{UInt8}) libpq
#@c Cint            PQgetCopyData (Ptr{PGconn}, Ptr{Ptr{UInt8}}, Cint) libpq

#### Results
@c Ptr{UInt8}   PQgetvalue (Ptr{PGresult}, Cint, Cint) libpq
@c Cint         PQgetisnull (Ptr{PGresult}, Cint, Cint) libpq
@c Void         PQclear (Ptr{PGresult},) libpq
# result fields
@c Cint         PQntuples (Ptr{PGresult},) libpq
@c Cint         PQnfields (Ptr{PGresult},) libpq
@c Cint         PQbinaryTuples (Ptr{PGresult},) libpq
@c Ptr{UInt8}   PQfname (Ptr{PGresult}, Cint) libpq
@c Cint         PQfnumber (Ptr{PGresult}, Ptr{UInt8}) libpq
@c Oid          PQftable (Ptr{PGresult}, Cint) libpq
@c Cint         PQftablecol (Ptr{PGresult}, Cint) libpq
@c Oid          PQftype (Ptr{PGresult}, Cint) libpq
#for update insert etc...
@c Ptr{UInt8}   PQcmdStatus (Ptr{PGresult},) libpq
@c Ptr{UInt8}   PQcmdTuples (Ptr{PGresult},) libpq

#### Escaping
@c Void         PQfreemem (Ptr{Void},) libpq
@c Ptr{UInt8}   PQescapeLiteral (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
#@c Ptr{UInt8} PQescapeIdentifier (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq

#### Canceling
@c Ptr{PGcancel} PQgetCancel (Ptr{PGconn},) libpq
@c Ptr{Void}     PQfreeCancel (Ptr{PGcancel}, ) libpq
@c Cint          PQcancel (Ptr{PGcancel}, Ptr{UInt8}, Cint) libpq

#### Misc
@c Cint PQprotocolVersion (Ptr{PGconn},) libpq
@c Cint PQserverVersion (Ptr{PGconn},) libpq
@c Cint PQlibVersion (Ptr{Void},) libpq
#@c Void PQreset (Ptr{PGconn},) libpq
#@c PGTransactionStatusType PQtransactionStatus (Ptr{PGconn},) libpq
#@c Ptr{Cuchar} PQescapeByteaConn (Ptr{PGconn}, Ptr{Cuchar}, Cint, Ptr{Cint}) libpq
#@c Ptr{Cuchar} PQunescapeBytea (Ptr{Cuchar}, Ptr{Cint}) libpq
#@c Ptr{Cuchar} PQescapeBytea (Ptr{Cuchar}, Cint, Ptr{Cint}) libpq

# return query as vector of vectors with all data types as strings
function bootstrap_query(ptr::Ptr{PGconn}, query::AbstractString)

    result = PQexec(ptr, query)
    ncols = PQnfields(result)
    nrows = PQntuples(result)
    a = Vector{Vector{UTF8String}}(nrows)
    for x in 1:nrows
        a[x] = [utf8(PQgetvalue(result, x-1, y-1)) for y in 1:ncols]
    end
    PQclear(result)
    return a
end

# not sure how to test
function interuptable_exec(ptr::Ptr{PGconn}, query::AbstractString)
    try
        res = PQexec(ptr, query)
    catch InteruptException
        cancel = PQgetCancel(ptr)
        msg = Array(UInt8, 256)
        status = PQcancel(cancel, msg, sizeof(msg))
        if status != 1
            error("cancel failed: $(bytestring(msg))")
        end
        PQfreeCancel(cancel)
        info("canceling statement due to user request")
        nothing
    end
end

ext: jl

max_stars_repo_name: NCarson/Postgres

max_issues_repo_name: NCarson/Postgres

max_forks_repo_name: NCarson/Postgres


text:
using Revise # reduce recompilation time
using Plots
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using StaticArrays
using Polyester
using MuladdMacro
using Setfield

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
include("../EntropyStableEuler.jl/src/logmean.jl")
using EntropyStableEuler
using EntropyStableEuler.Fluxes2D

@muladd begin

const TOL     = 1e-15
const POSTOL  = 1e-14   # Tolerance of positivity
const Nc      = 4       # number of components
const WALLPT  = 1.0/6.0 # Starting point of the wall in x direction
# const XLENGTH = 4.0     # length of domain TODO: hardcoded now

"Approximation parameters"
const N   = 3
const K1D = 200
const T   = 0.1

const Nq  = N
const Np  = (N+1)*(N+1)
const Nfp = 4*(N+1)

const NUM_THREADS = Threads.nthreads()

# Initial condition 2D shocktube
const γ = 1.4
const rhoL  = 8.0
const uL    = 8.25*cos(pi/6)
const vL    = -8.25*sin(pi/6)
const pL    = 116.5
const rhouL = rhoL*uL
const rhovL = rhoL*vL
const EL    = pL/(γ-1)+.5*rhoL*(uL^2+vL^2)
const betaL = rhoL/(2*pL)
const rhoR  = 1.4
const uR    = 0.0
const vR    = 0.0
const pR    = 1.0
const rhouR = rhoR*uR
const rhovR = rhoR*vR
const ER    = pR/(γ-1)+.5*rhoR*(uR^2+vR^2)
const betaR = rhoR/(2*pR)
const SHOCKSPD = 10.0/cos(pi/6)

"Mesh related variables"
VX, VY, EToV = uniform_quad_mesh(4*K1D,K1D)
@. VX = (VX+1)*2
@. VY = (VY+1)/2

"Initialize reference element"
rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
@unpack Vf,Dr,Ds,LIFT,Nfaces = rd

md = init_mesh((VX,VY),EToV,rd)
@unpack sJ,xf,yf,mapM,mapP,mapB,nxJ,nyJ,x,y = md
# TODO: assume uniform mesh, hardcoded geometric factors
const K   = 4*K1D*K1D
const J   = 1.0/K
const Jf  = 1.0/K1D/2
const rxJ = 2*K1D*J
const syJ = 2*K1D*J

xb,yb = (x->x[mapB]).((xf,yf))
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

rq,wq = gauss_lobatto_quad(0,0,Nq)
M = diagm(wq)
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# 1D operators
M1D = Mlump
M1D_inv = Mlump_inv
Q1D = M*Dr
B1D = zeros(N+1,N+1)
B1D[1,1] = -1
B1D[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B1D*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Q01D = S0+1/2*B1D

# Drop zeros
Q1D = Matrix(droptol!(sparse(Q1D),TOL))
Q01D = Matrix(droptol!(sparse(Q01D),TOL))
B1D = Matrix(droptol!(sparse(B1D),TOL))
S01D = Matrix(droptol!(sparse(S0),TOL))
S1D = Matrix(droptol!(sparse((Q1D-Q1D')/2),TOL))

# Tensor product operators
M = droptol!(sparse(kron(M1D,M1D)),TOL)
Minv = spdiagm(0 => 1 ./ diag(M))
Qr = droptol!(sparse(kron(M1D,Q1D)),TOL)
Qs = droptol!(sparse(kron(Q1D,M1D)),TOL)
Qr0 = droptol!(sparse(kron(M1D,Q01D)),TOL)
Qs0 = droptol!(sparse(kron(Q01D,M1D)),TOL)
Sr = droptol!(sparse((Qr-Qr')/2),TOL)
Ss = droptol!(sparse((Qs-Qs')/2),TOL)
S0r = droptol!(sparse(kron(M1D,S01D)),TOL)
S0s = droptol!(sparse(kron(S01D,M1D)),TOL)
Br = droptol!(sparse(Qr+Qr'),TOL)
Bs = droptol!(sparse(Qs+Qs'),TOL)

# Convert diagonal matrices to vectors
Br   = Array(diag(Br))
Bs   = Array(diag(Bs))
M    = Array(diag(M))
Minv = Array(diag(Minv))

Fmask  = [1:N+1; (N+1):(N+1):Np; Np:-1:Np-N; Np-N:-(N+1):1]
Fxmask = [(N+2):(2*N+2); (3*N+4):(4*N+4)]
Fymask = [1:(N+1); (2*N+3):(3*N+3)]

# Boundary conditions
inflow   = mapB[findall(@. (abs(xb) < TOL) | ((xb < 1/6) & (abs(yb) < TOL)))]
outflow  = mapB[findall(@. abs(xb-4.) < TOL)]
topflow  = mapB[findall(@. abs(yb-1.) < TOL)]
wall     = mapB[findall(@. (xb >= 1/6) & (abs(yb) < TOL))]
nx_wall  = nxJ[wall]./sJ[wall]
ny_wall  = nyJ[wall]./sJ[wall]
const TOP_INIT = (1+sqrt(3)/6)/sqrt(3)

@inline function get_valP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)
    inflow  = ((abs(xM) < TOL) & (abs(yM-1.) > TOL)) | ((xM <= WALLPT) & (abs(yM) < TOL))
    outflow = ((abs(xM-4.0) < TOL) & (abs(yM) > TOL) & (abs(yM-1.0) > TOL))
    topflow = (abs(yM-1.0) < TOL)
    wall    = ((xM > WALLPT) & (abs(yM) < TOL))
    has_bc  = (inflow | outflow | topflow | wall)

    if inflow
        rhoP   = rhoL
        rhouP  = rhouL
        rhovP  = rhovL
        EP     = EL
        fx_1_P = rhoL*uL
        fx_2_P = rhoL*uL^2+pL
        fx_3_P = rhoL*uL*vL
        fx_4_P = uL*(EL+pL)
        fy_1_P = rhoL*vL
        fy_2_P = rhoL*uL*vL
        fy_3_P = rhoL*vL^2+pL
        fy_4_P = vL*(EL+pL)

    elseif outflow
        rhoP   = U[1,iM,k]
        rhouP  = U[2,iM,k]
        rhovP  = U[3,iM,k]
        EP     = U[4,iM,k]
        fx_1_P = f_x[1,iM,k]
        fx_2_P = f_x[2,iM,k]
        fx_3_P = f_x[3,iM,k]
        fx_4_P = f_x[4,iM,k]
        fy_1_P = f_y[1,iM,k]
        fy_2_P = f_y[2,iM,k]
        fy_3_P = f_y[3,iM,k]
        fy_4_P = f_y[4,iM,k]

    elseif wall
        # TODO: we assume the normals are [0;-1] here
        # Un = -vM
        # Ut = -uM
        rhoP   = U[1,iM,k]
        rhouP  = rhoP*uM   # rhoP = rhoM, TODO: refactor
        rhovP  = -rhoP*vM
        uP     = rhouP/rhoP
        vP     = rhovP/rhoP
        EP     = U[4,iM,k]
        pP     = pfun(rhoP,rhouP,rhovP,EP)
        fx_1_P = rhoP*uP
        fx_2_P = rhoP*uP^2+pP
        fx_3_P = rhoP*uP*vP
        fx_4_P = uP*(EP+pP)
        fy_1_P = rhoP*vP
        fy_2_P = rhoP*uP*vP
        fy_3_P = rhoP*vP^2+pP
        fy_4_P = vP*(EP+pP)

    elseif topflow
        breakpoint = TOP_INIT+t*SHOCKSPD
        if xM < breakpoint
            rhoP   = rhoL
            rhouP  = rhouL
            rhovP  = rhovL
            EP     = EL
            fx_1_P = rhoL*uL
            fx_2_P = rhoL*uL^2+pL
            fx_3_P = rhoL*uL*vL
            fx_4_P = uL*(EL+pL)
            fy_1_P = rhoL*vL
            fy_2_P = rhoL*uL*vL
            fy_3_P = rhoL*vL^2+pL
            fy_4_P = vL*(EL+pL)
        else
            rhoP   = rhoR
            rhouP  = rhouR
            rhovP  = rhovR
            EP     = ER
            fx_1_P = rhoR*uR
            fx_2_P = rhoR*uR^2+pR
            fx_3_P = rhoR*uR*vR
            fx_4_P = uR*(ER+pR)
            fy_1_P = rhoR*vR
            fy_2_P = rhoR*uR*vR
            fy_3_P = rhoR*vR^2+pR
            fy_4_P = vR*(ER+pR)
        end
    else                         # if not on the physical boundary
        gP = mapP[i,k]           # exterior global face node number
        kP = fld1(gP,Nfp)        # exterior element number
        iP = Fmask[mod1(gP,Nfp)] # exterior node number

        rhoP  = U[1,iP,kP]
        rhouP = U[2,iP,kP]
        rhovP = U[3,iP,kP]
        EP    = U[4,iP,kP]
        fx_1_P = f_x[1,iP,kP]
        fx_2_P = f_x[2,iP,kP]
        fx_3_P = f_x[3,iP,kP]
        fx_4_P = f_x[4,iP,kP]
        fy_1_P = f_y[1,iP,kP]
        fy_2_P = f_y[2,iP,kP]
        fy_3_P = f_y[3,iP,kP]
        fy_4_P = f_y[4,iP,kP]
    end

    return rhoP,rhouP,rhovP,EP,fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P,has_bc
end

# initial condition
at_left(x,y) = y-sqrt(3)*x+sqrt(3)/6 > 0.0
U = zeros(Float64,Nc,Np,K)
for k = 1:K
    for i = 1:Np
        if at_left(x[i,k],y[i,k])
            U[1,i,k] = rhoL
            U[2,i,k] = rhouL
            U[3,i,k] = rhovL
            U[4,i,k] = EL
        else
            U[1,i,k] = rhoR
            U[2,i,k] = rhouR
            U[3,i,k] = rhovR
            U[4,i,k] = ER
        end
    end
end

@inline function pfun(rho,rhou,E)
    return (γ-1)*(E-.5*rhou^2/rho)
end

@inline function pfun(rho,rhou,rhov,E)
    return (γ-1)*(E-.5*(rhou^2+rhov^2)/rho)
end

@inline function wavespeed_1D(rho,rhou,E)
    p = pfun(rho,rhou,E)
    return abs(rhou/rho) + sqrt(γ*p/rho)
end

@inline function logmean(aL,aR,logL,logR)

    # "from: Entropy stable num. approx. for the isothermal and polytropic Euler"

    da = aR-aL;
    aavg = .5*(aR+aL);
    f = da/aavg;
    v = f^2;
    if abs(f)<1e-4
        # numbers assume the specific value γ = 1.4
        return aavg*(1 + v*(-.2-v*(.0512 - v*0.026038857142857)))
    else
        return -da/(logL-logR)
    end
end

@inline function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                                 rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    FyS1 = rholog*vavg
    FyS2 = FxS3
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FxS1,FxS2,FxS3,FxS4,FyS1,FyS2,FyS3,FyS4
end

@inline function limiting_param(rhoL,rhouL,rhovL,EL,rhoP,rhouP,rhovP,EP)
    # L - low order, P - P_ij

    l = 1.0
    # Limit density
    if rhoL + rhoP < -TOL
        l = max((-rhoL+POSTOL)/rhoP, 0)
    end

    # limiting internal energy (via quadratic function)
    a = rhoP*EP-(rhouP^2+rhovP^2)/2.0
    b = rhoP*EL+rhoL*EP-rhouL*rhouP-rhovL*rhovP
    c = rhoL*EL-(rhouL^2+rhovL^2)/2.0-POSTOL

    l_eps_ij = 1.0
    if b^2-4*a*c >= 0
        r1 = (-b+sqrt(b^2-4*a*c))/(2*a)
        r2 = (-b-sqrt(b^2-4*a*c))/(2*a)
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

# Preallocation
f_x     = zeros(Float64,size(U))
f_y     = zeros(Float64,size(U))
rholog  = zeros(Float64,Np,K)
betalog = zeros(Float64,Np,K)
U_low   = zeros(Float64,Nc,Np,NUM_THREADS)
F_low   = zeros(Float64,Nc,Np,Np,NUM_THREADS)
F_high  = zeros(Float64,Nc,Np,Np,NUM_THREADS)
F_P     = zeros(Float64,Nc,Nfp,NUM_THREADS)
L       =  ones(Float64,Np,Np,NUM_THREADS)

const DEBUGELEM = 3
function rhs_IDP_fixdt!(U,rhsU,t,prealloc,ops,geom)
    f_x,f_y,U_low,F_low,F_high,F_P,L,rholog,betalog = prealloc
    S0r,S0s,Sr,Ss,Br_halved,Bs_halved,MJ_inv,coeff_arr = ops
    mapP,Fmask,Fxmask,Fymask,x,y = geom

    # Preallocate inviscid flux TODO: shall we?
    # Current time: 107 ms
    # Current time: 27.5 ms
    # TODO: rewrite arithmetic expressions
    # TODO: BC
    # TODO: BC dissipation
    # TODO: @inbounds @turbo
    fill!(rhsU, 0.0)   # TODO: fill rhsU with zeros?
    for e = 1:K
        for i = 1:Np
            rho  = U[1,i,e]
            rhou = U[2,i,e]
            rhov = U[3,i,e]
            E    = U[4,i,e]
            p          = pfun(rho,rhou,rhov,E)
            f_x[1,i,e] = rhou
            f_x[2,i,e] = rhou^2/rho+p
            f_x[3,i,e] = rhou*rhov/rho
            f_x[4,i,e] = E*rhou/rho+p*rhou/rho
            f_y[1,i,e] = rhov
            f_y[2,i,e] = rhou*rhov/rho
            f_y[3,i,e] = rhov^2/rho+p
            f_y[4,i,e] = E*rhov/rho+p*rhov/rho
            rholog[i,e]  = log(rho)
            betalog[i,e] = log(rho/(2*p))
        end
    end

    # =====================
    # Loop through elements
    # =====================
    for  k = 1:K
        tid = Threads.threadid()

        # ===========================================
        # Reinitialize arrays TODO: maybe unnecessary
        # ===========================================
        fill!(F_low , 0.0)
        fill!(F_high, 0.0)
        fill!(U_low , 0.0)
        fill!(F_P   , 0.0)
        fill!(L     , 1.0)

        # Calculate low order algebraic flux
        for j = 2:Np
            rho_j  = U[1,j,k]
            rhou_j = U[2,j,k]
            rhov_j = U[3,j,k]
            E_j    = U[4,j,k]
            fx_1_j = f_x[1,j,k]
            fx_2_j = f_x[2,j,k]
            fx_3_j = f_x[3,j,k]
            fx_4_j = f_x[4,j,k]
            fy_1_j = f_y[1,j,k]
            fy_2_j = f_y[2,j,k]
            fy_3_j = f_y[3,j,k]
            fy_4_j = f_y[4,j,k]
            for i = 1:j-1
                S0r_ij = S0r[i,j]
                S0s_ij = S0s[i,j]
                rho_i  = U[1,i,k]
                rhou_i = U[2,i,k]
                rhov_i = U[3,i,k]
                E_i    = U[4,i,k]
                fx_1_i = f_x[1,i,k]
                fx_2_i = f_x[2,i,k]
                fx_3_i = f_x[3,i,k]
                fx_4_i = f_x[4,i,k]
                fy_1_i = f_y[1,i,k]
                fy_2_i = f_y[2,i,k]
                fy_3_i = f_y[3,i,k]
                fy_4_i = f_y[4,i,k]

                if S0r_ij != 0 || S0s_ij != 0
                    n_ij_norm = sqrt(rxJ_sq*S0r_ij^2+syJ_sq*S0s_ij^2)
                    n_ij_x = rxJ*S0r_ij/n_ij_norm
                    n_ij_y = syJ*S0s_ij/n_ij_norm
                    λ_i = wavespeed_1D(rho_i,n_ij_x*rhou_i+n_ij_y*rhov_i,E_i)
                    λ_j = wavespeed_1D(rho_j,n_ij_x*rhou_j+n_ij_y*rhov_j,E_j)
                    λ_ij = max(λ_i,λ_j)*n_ij_norm

                    FL1 = (rxJ*S0r_ij*(fx_1_i+fx_1_j)
                          +syJ*S0s_ij*(fy_1_i+fy_1_j)
                          -λ_ij*(rho_j-rho_i) )
                    FL2 = (rxJ*S0r_ij*(fx_2_i+fx_2_j)
                          +syJ*S0s_ij*(fy_2_i+fy_2_j)
                          -λ_ij*(rhou_j-rhou_i) )
                    FL3 = (rxJ*S0r_ij*(fx_3_i+fx_3_j)
                          +syJ*S0s_ij*(fy_3_i+fy_3_j)
                          -λ_ij*(rhov_j-rhov_i) )
                    FL4 = (rxJ*S0r_ij*(fx_4_i+fx_4_j)
                          +syJ*S0s_ij*(fy_4_i+fy_4_j)
                          -λ_ij*(E_j-E_i) )

                    F_low[1,i,j,tid] = FL1
                    F_low[2,i,j,tid] = FL2
                    F_low[3,i,j,tid] = FL3
                    F_low[4,i,j,tid] = FL4

                    F_low[1,j,i,tid] = -FL1
                    F_low[2,j,i,tid] = -FL2
                    F_low[3,j,i,tid] = -FL3
                    F_low[4,j,i,tid] = -FL4
                end
            end
        end

        # Calculate high order algebraic flux
        for j = 2:Np
            rho_j     = U[1,j,k]
            u_j       = U[2,j,k]/rho_j
            v_j       = U[3,j,k]/rho_j
            beta_j    = rho_j/(2*pfun(rho_j,rho_j*u_j,rho_j*v_j,U[4,j,k]))
            rholog_j  = rholog[j,k]
            betalog_j = betalog[j,k]
            for i = 1:j-1
                Sr_ij     = Sr[i,j]
                Ss_ij     = Ss[i,j]
                rho_i     = U[1,i,k]
                u_i       = U[2,i,k]/rho_i
                v_i       = U[3,i,k]/rho_i
                beta_i    = rho_i/(2*pfun(rho_i,rho_i*u_i,rho_i*v_i,U[4,i,k]))
                rholog_i  = rholog[i,k]
                betalog_i = betalog[i,k]
                if Sr_ij != 0.0 || Ss_ij != 0.0
                    Fx1,Fx2,Fx3,Fx4,Fy1,Fy2,Fy3,Fy4 = euler_fluxes_2D(rho_i,u_i,v_i,beta_i,rholog_i,betalog_i,
                                                                      rho_j,u_j,v_j,beta_j,rholog_j,betalog_j)

                    FH1 = rxJ_db*Sr_ij*Fx1+syJ_db*Ss_ij*Fy1
                    FH2 = rxJ_db*Sr_ij*Fx2+syJ_db*Ss_ij*Fy2
                    FH3 = rxJ_db*Sr_ij*Fx3+syJ_db*Ss_ij*Fy3
                    FH4 = rxJ_db*Sr_ij*Fx4+syJ_db*Ss_ij*Fy4

                    F_high[1,i,j,tid] = FH1
                    F_high[2,i,j,tid] = FH2
                    F_high[3,i,j,tid] = FH3
                    F_high[4,i,j,tid] = FH4

                    F_high[1,j,i,tid] = -FH1
                    F_high[2,j,i,tid] = -FH2
                    F_high[3,j,i,tid] = -FH3
                    F_high[4,j,i,tid] = -FH4
                end
            end
        end

        # Calculate interface fluxes
        for i = 1:Nfp
            Br_ii_halved = Br_halved[Fmask[i]]
            Bs_ii_halved = Bs_halved[Fmask[i]]

            iM    = Fmask[i]
            xM    = x[iM,k]
            yM    = y[iM,k]
            rhoM  = U[1,iM,k]
            rhouM = U[2,iM,k]
            rhovM = U[3,iM,k]
            EM    = U[4,iM,k]
            uM    = rhouM/rhoM
            vM    = rhovM/rhoM
            fx_1_M = f_x[1,iM,k]
            fx_2_M = f_x[2,iM,k]
            fx_3_M = f_x[3,iM,k]
            fx_4_M = f_x[4,iM,k]
            fy_1_M = f_y[1,iM,k]
            fy_2_M = f_y[2,iM,k]
            fy_3_M = f_y[3,iM,k]
            fy_4_M = f_y[4,iM,k]

            rhoP,rhouP,rhovP,EP,fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P,has_bc = get_valP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)

            # flux in x direction
            if i in Fxmask
                λM = wavespeed_1D(rhoM,rhouM,EM)
                λP = wavespeed_1D(rhoP,rhouP,EP)
                if has_bc
                    λ = 0
                else
                    λ  = max(λM,λP)*abs(Br_ii_halved)
                end

                F_P[1,i,tid] = (Jf*Br_ii_halved*(fx_1_M+fx_1_P)
                               -λ*(rhoP-rhoM) )
                F_P[2,i,tid] = (Jf*Br_ii_halved*(fx_2_M+fx_2_P)
                               -λ*(rhouP-rhouM) )
                F_P[3,i,tid] = (Jf*Br_ii_halved*(fx_3_M+fx_3_P)
                               -λ*(rhovP-rhovM) )
                F_P[4,i,tid] = (Jf*Br_ii_halved*(fx_4_M+fx_4_P)
                               -λ*(EP-EM) )
            end

            # flux in y direction
            if i in Fymask
                λM = wavespeed_1D(rhoM,rhovM,EM)
                λP = wavespeed_1D(rhoP,rhovP,EP)
                if has_bc
                    λ = 0
                else
                    λ  = max(λM,λP)*abs(Bs_ii_halved)
                end

                F_P[1,i,tid] = (Jf*Bs_ii_halved*(fy_1_M+fy_1_P)
                               -λ*(rhoP-rhoM) )
                F_P[2,i,tid] = (Jf*Bs_ii_halved*(fy_2_M+fy_2_P)
                               -λ*(rhouP-rhouM) )
                F_P[3,i,tid] = (Jf*Bs_ii_halved*(fy_3_M+fy_3_P)
                               -λ*(rhovP-rhovM) )
                F_P[4,i,tid] = (Jf*Bs_ii_halved*(fy_4_M+fy_4_P)
                               -λ*(EP-EM) )
            end
        end

        # Calculate low order solution
        for i = 1:Np
            for j = 1:Np
                for c = 1:Nc
                    U_low[c,i,tid] = U_low[c,i,tid] + F_low[c,i,j,tid]
                end
            end
        end

        for i = 1:Nfp
            iM = Fmask[i]
            for c = 1:Nc
                U_low[c,iM,tid] = U_low[c,iM,tid] + F_P[c,i,tid]
            end
        end

        for c = 1:Nc
            for i = 1:Np
                U_low[c,i,tid] = U[c,i,k] - dt*MJ_inv[i]*U_low[c,i,tid]
            end
        end

        # Calculate limiting parameters
        for i = 1:Np
            coeff = coeff_arr[i]
            rhoL  = U_low[1,i,tid]
            rhouL = U_low[2,i,tid]
            rhovL = U_low[3,i,tid]
            EL    = U_low[4,i,tid]
            for j = 1:Np
                if i != j
                    rhoP  = coeff*(F_low[1,i,j,tid]-F_high[1,i,j,tid])
                    rhouP = coeff*(F_low[2,i,j,tid]-F_high[2,i,j,tid])
                    rhovP = coeff*(F_low[3,i,j,tid]-F_high[3,i,j,tid])
                    EP    = coeff*(F_low[4,i,j,tid]-F_high[4,i,j,tid])
                    L[i,j,tid] = limiting_param(rhoL,rhouL,rhovL,EL,rhoP,rhouP,rhovP,EP)
                end
            end
        end

        # Elementwise limiting
        l_e = 1.0
        for j = 1:Np
            for i = 1:Np
                if i != j
                    l_e = min(l_e,L[i,j,tid])
                end
            end
        end
        for j = 1:Np
            for i = 1:Np
                # TODO: debug
                L[i,j,tid] = 0.0#l_e
            end
        end

        # Accumulate RHS
        # for j = 1:Np
        #     for i = 1:Np
        #         for c = 1:Nc
        #             rhsU[c,i,k] = rhsU[c,i,k] + (L[i,j,tid]-1)*F_low[c,i,j,tid]-L[i,j,tid]*F_high[c,i,j,tid]
        #         end
        #     end
        # end
        # for i = 1:Nfp
        #     for c = 1:Nc
        #         rhsU[c,Fmask[i],k] = rhsU[c,Fmask[i],k] - F_P[c,i,tid]
        #     end
        # end

        for c = 1:Nc
            for i = 1:Np
                for j = 1:Np
                    rhsU[c,i,k] = rhsU[c,i,k] + (L[i,j,tid]-1)*F_low[c,i,j,tid]-L[i,j,tid]*F_high[c,i,j,tid]
                    # rhsU[c,i,k] = rhsU[c,i,k] - F_low[c,i,j,tid]
                end
            end
            for i = 1:Nfp
                rhsU[c,Fmask[i],k] = rhsU[c,Fmask[i],k] - F_P[c,i,tid]
            end
        end

        # if k == DEBUGELEM
        #     for c = 1:4
        #         println(c)
        #         println("=== F_low ===", c)
        #         display(F_low[c,:,:,1])
        #         println("=== F_P ===", c)
        #         display(F_P[c,:,1])
        #         println("=== U_low ===", c)
        #         display(U_low[c,:,1])
        #         println("=== rhsU ===", c)
        #         display(rhsU[c,:,DEBUGELEM])
        #     end
        # end
    end

    for k = 1:K
        for i = 1:Np
            for c = 1:Nc
                rhsU[c,i,k] = MJ_inv[i]*rhsU[c,i,k]
            end
        end
    end
end

t = 0.0
const dt = 1e-4
rhsU = zeros(size(U))

# Define some constants for optimizations
const rxJ_sq = rxJ^2
const syJ_sq = syJ^2
const rxJ_db = 2*rxJ
const syJ_db = 2*syJ
MJ_inv    = Minv./J
# Br_halved = Br./2.0
# Bs_halved = Bs./2.0
Br_halved = -sum(S0r,dims=2)
Bs_halved = -sum(S0s,dims=2)
coeff_arr = dt*(Np-1).*Minv*J

prealloc = (f_x,f_y,U_low,F_low,F_high,F_P,L,rholog,betalog)
ops      = (S0r,S0s,Sr,Ss,Br_halved,Bs_halved,MJ_inv,coeff_arr)
geom     = (mapP,Fmask,Fxmask,Fymask,x,y)

t = 0.0
#@btime rhs_IDP_fixdt!($U,$rhsU,$t,$prealloc,$ops,$geom)
@profiler rhs_IDP_fixdt!(U,rhsU,t,prealloc,ops,geom)
#rhs_IDP_fixdt!(U,rhsU,t,prealloc,ops,geom)

# for c = 1:4
#     println("=== rhsU ===", c)
#     display(rhsU[c,:,:])
# end

# # Time stepping
# i = 1
# t = 0.0
# rhsU = zeros(size(U))
# resW = zeros(size(U))

# # while t < T
# while i < 2
#     # TODO: fix
#     # dt = min(1e-4,T-t)

#     # println("===== First step =====")
#     # for c = 1:4
#     #     println("=== U ===", c)
#     #     display(U[c,:,:])
#     # end
#     rhs_IDP_fixdt!(U,rhsU,t,prealloc,ops,geom)
#     @. resW = U + dt*rhsU
#     println("===== Second step =====")
#     for c = 1:4
#         println("=== resW ===", c)
#         display(resW[c,:,:])
#     end
#     rhs_IDP_fixdt!(resW,rhsU,t,prealloc,ops,geom)
#     println("===== Second RHSU =====")
#     for c = 1:4
#         println("=== rhsU ===", c)
#         display(rhsU[c,:,:])
#     end
#     # println("===== Look at first resW =====")
#     # for c = 1:4
#     #     println("=== resW ===", c)
#     #     display(resW[c,:,:])
#     # end
#     @. resW = resW+dt*rhsU
#     # println("===== Look at second resW =====")
#     # for c = 1:4
#     #     println("=== resW ===", c)
#     #     display(resW[c,:,:])
#     # end
#     @. resW = 0.75*U+0.25*resW
#     # println("===== Look at third resW =====")
#     # for c = 1:4
#     #     println("=== resW ===", c)
#     #     display(resW[c,:,:])
#     # end
#     # println("===== Third step =====")
#     # for c = 1:4
#     #     println("=== resW ===", c)
#     #     display(resW[c,:,:])
#     # end
#     rhs_IDP_fixdt!(resW,rhsU,t,prealloc,ops,geom)
#     @. resW = resW+dt*rhsU
#     @. U = 1/3*U + 2/3*resW
#     # println("===== Last step =====")
#     # for c = 1:4
#     #     println("=== U ===", c)
#     #     display(U[c,:,:])
#     # end
#     global t = t + dt
#     println("Current time $t with time step size $dt, and final time $T, at step $i")
#     global i = i + 1
# end

# # rp,sp = equi_nodes_2D(10)
# # @unpack VDM = rd
# # Vp = vandermonde_2D(N,rp,sp)/VDM
# # gr(aspect_ratio=:equal,legend=false,markerstrokewidth=0,markersize=2)
# # xp = Vp*x
# # yp = Vp*y
# # rho = Vp*U[1,:,:]

# # # scatter(xp,yp,zcolor=rho,camera=(0,90))
# # # savefig("~/Desktop/tmp.png")

end # muladd

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
using Interpolations, ColorVectorSpace, ImageFiltering, Images, ScatteredInterpolation

#TODO edit doc
"""
    warp_img(img::Image, dx, dy; border_strat::Symbol=:replicate)::Image
    warp_img(img::Image, dx, dy, filling_img::Image)::Image

Warp the image `img` by `dx` in *x* direction and by `dy` in *y* direction.
`border_strat` indicates how to act when the resulting coordinate outside of bounds of `img`.

# Possible variations:
- `border_strat == :replicate`: fill with the closest pixel value.
- `border_strat == :zeros`: fill with zeros.
- supply a `filling_img` as 3rd argument: fill with the provided image `filling_img`

See also: [`showflow`](@ref), [`imgshowflow`](@ref), [`Flow`](@ref)
"""
function warp_img(img, dx, dy; border_strat::Symbol=:replicate)
    if border_strat == :replicate
        return img_warp_replicate(img, dx, dy)
    elseif border_strat == :zeros
        return warp_img_zeros(img, dx, dy)
    end
end

function warp_img(img::Image, dx, dy, filling_img::Image)::Image
    return warp_img_with_filling(img, dx, dy, filling_img)
end

function warp_img_zeros(img::Image, dx, dy)::Image
    itp = Interpolations.interpolate(img, BSpline(Linear()))
    etp = Interpolations.extrapolate(itp, 0)
    inds = indices_spatial(img)
    rng = extrema.(inds)
    imgw = similar(img, eltype(itp))
    for I in CartesianIndices(inds)
        y, x = Tuple(I) .+ (dy[I], dx[I])
        imgw[I] = etp(y, x)
    end
    return imgw
end

function warp_img_with_filling(img::Image, dx, dy, filling_img::Image)::Image
    itp = Interpolations.interpolate(img, BSpline(Linear()))
    inds = indices_spatial(img)
    rng = extrema.(inds)
    imgw = similar(img, eltype(itp))
    for I in CartesianIndices(inds)
        imgw[I] = interplate_or_fill(I, dx[I], dy[I], rng, itp, filling_img)
    end
    return imgw
end

@inline function interplate_or_fill(ind, dxi, dyi, rng, itp, filling_img)
    y, x = (ind[1]+dyi, ind[2]+dxi)
    if is_in_bounds(y, rng[1]...) && is_in_bounds(x, rng[2]...)
        return itp(y, x)
    else
        return filling_img[ind]
    end
end

@inline function is_in_bounds(num, lo, hi)
    return lo <= num <= hi
end

function img_warp_replicate(img, dx, dy)
    itp = Interpolations.interpolate(img, BSpline(Linear()))
    inds = indices_spatial(img)
    rng = extrema.(inds)
    imw = similar(img, eltype(itp))
    for I in CartesianIndices(inds)
        dxi, dyi = dx[I], dy[I]
        y, x = clamp(I[1]+dyi, rng[1]...), clamp(I[2]+dxi, rng[2]...)
        imw[I] = itp(y, x)
    end
    return imw
end

meshgrid(x, y) = [repeat(x, outer=length(y)) repeat(y, inner=length(x))]
meshgrid(x::Real, y::Real) = [repeat(1:x, outer=y) repeat(1:y, inner=x)]

"""
interpolate_flow(flow_at_inds,
                 inds::Array{CartesianIndex{2}, 1},
                 flow_size;
                 method::Symbol=:quad,
                 kwargs=Dict())

Interpolate a complex dispalcment field of size `flow_size` from an array of displacement vectors `flow_at_inds`, that are at the locations `inds`.
You can choose between two methods `method=:quad` which fits the displacement vectors into a global quadratic displacement field
and `method=:rbf` which uses an rbf model. Additional keyword arguments of the chosen method can be passed as `kwargs`.

See also: [`Flow`](@ref), [`interpolate_flow_quad`](@ref), [`interpolate_flow_rbf`](@ref)
"""
function interpolate_flow(flow_at_inds,
                          inds::Array{CartesianIndex{2}, 1},
                          flow_size;
                          method::Symbol=:quad,
                          kwargs=Dict())

if method == :quad
    generated_flow = interpolate_flow_quad(flow_at_inds, inds, flow_size)
elseif method == :rbf
    generated_flow = interpolate_flow_rbf(flow_at_inds, inds, flow_size; kwargs...)
end
    return generated_flow
end

#TODO edit docs
"""
    interpolate_flow_rbf(flow_at_inds,
                         inds::Array{CartesianIndex{2}, 1},
                         flow_size,
                         rbf_method::T=Multiquadratic(2)) where {T <: ScatteredInterpolation.RadialBasisFunction}

Interpolate a complex displacement field of size `flow_size` from displacement vectors `flow_at_inds` at `inds`
using the a Multiquadratic RBF model with the parameter `ε` and using only the values at `inds` for the interpolation.

# Examples of `method` from ScatteredInterpolation:
Multiquadratic(ɛ = 1)
```math
ϕ(r) = \\sqrt{1 + (ɛr)^2}
```
Polyharmonic(k = 1)
```math
ϕ(r) = r^k, k = 1, 3, 5, ...
\\\\
ϕ(r) = r^k ln(r), k = 2, 4, 6, ...
```
(See ScatteredInterpolation for more.)

See also: [`showflow`](@ref), [`Flow`](@ref), [`interpolate_flow`](@ref), [`interpolate_flow_quad`](@ref)
"""
function interpolate_flow_rbf(flow_at_inds,
                              inds::Array{CartesianIndex{2}, 1},
                              flow_size,
                              rbf_method::T=Multiquadratic(2)) where {T <: ScatteredInterpolation.RadialBasisFunction}

    non_nan_locs = filter(x -> !isnan(flow_at_inds[x]), 1:length(flow_at_inds))
    @assert length(non_nan_locs) != 0

    # interpolate
    itp = ScatteredInterpolation.interpolate(rbf_method, inds_to_points(inds[non_nan_locs]), flow_at_inds[non_nan_locs]);
    gridPoints = meshgrid(flow_size...)'
    interpolated = ScatteredInterpolation.evaluate(itp, gridPoints)

    gridded = reshape(interpolated, flow_size)

    return gridded
end

"""
    interpolate_flow_quad(flow_at_inds, inds, flow_size)

Interpolate a complex displacement field of size `flow_size` from displacement vectors `flow_at_inds` at `inds`
using a global quadratic model.

The fit to this model is made by minimizing the following:

```math
\\min_{a_{k}} \\sum_{x, y ∈ \\mathrm{inds}}\\left|u_{\\mathrm{quad}}(x, y)-u_{\\mathrm{flow\\_at\\_inds}}(x, y)\\right|^{2}
```,

where:

```math
u_{\\mathrm{quad}}(x, y) = a_1 + a_2x + a_3y + a_4x^2 + a_5y^2 + a_6xy
```

See also: [`showflow`](@ref), [`Flow`](@ref), [`interpolate_flow`](@ref), [`interpolate_flow_rbf`](@ref)
"""
function interpolate_flow_quad(flow_at_inds, inds, flow_size)
    A = make_A(inds)
    coeffs = qr(A, Val(true))\flow_at_inds
    estim_flow = make_flow2(coeffs, flow_size)
end

function make_flow2(coeffs, flow_size)
    # flow = similar(coeffs, ComplexF64, flow_size)
    flow = Array{ComplexF64}(undef, flow_size)
    xs = collect(1:flow_size[1])
    if flow_size[1] == flow_size[2]
        ys = xs
    else
        ys = collect(1:flow_size[2])
    end
    flow .= coeffs[1]
    flow .+= (xs .* coeffs[2])
    transpose(flow) .+= (ys .* coeffs[3])
    flow .+= ((xs.^2) .* coeffs[4])
    transpose(flow) .+= ((ys.^2) .* coeffs[5])
    flow .+= ((repeat(xs, 1, flow_size[2]) .* ys') .* coeffs[6])
    return flow
end

function make_A(inds)
    points = transpose(inds_to_points(inds))
    A = hcat(ones(length(inds)), points[:, 1], points[:,2], points[:,1].^2, points[:,2].^2, points[:,1] .* points[:,2])
end


text:
parent_type(::Type{AlgAssAbsOrdElem{S, T}}) where {S, T} = AlgAssAbsOrd{S, T}

parent_type(::AlgAssAbsOrdElem{S, T}) where {S, T} = AlgAssAbsOrd{S, T}

@inline parent(x::AlgAssAbsOrdElem) = x.parent

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::AlgAssAbsOrdElem{S, T}, y::AlgAssAbsOrdElem{S, T}) where {S, T}
  return parent(x) === parent(y)
end

################################################################################
#
#  Parent object overloading
#
################################################################################

(O::AlgAssAbsOrd{S, T})(a::T, check::Bool = true) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    !x && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), y)
  else
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a))
  end
end

(O::AlgAssAbsOrd{S, T})(a::T, arr::Vector{fmpz}, check::Bool = false) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    (!x || arr != y) && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), y)
  else
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), deepcopy(arr))
  end
end

(O::AlgAssAbsOrd{S, T})(arr::Vector{fmpz}) where {S, T} = begin
  M = basis_mat(O, copy = false)
  N = matrix(FlintZZ, 1, degree(O), arr)
  NM = N*M
  x = elem_from_mat_row(algebra(O), NM.num, 1, NM.den)
  return AlgAssAbsOrdElem{S, T}(O, x, deepcopy(arr))
end

(O::AlgAssAbsOrd{S, T})(a::AlgAssAbsOrdElem{S, T}, check::Bool = true) where {S, T} = begin
  b = elem_in_algebra(a) # already a copy
  if check
    (x, y) = _check_elem_in_order(b, O)
    !x && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, b, y)
  else
    return AlgAssAbsOrdElem{S, T}(O, b)
  end
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdElem, dict::IdDict)
  b = parent(a)()
  b.elem_in_algebra = Base.deepcopy_internal(a.elem_in_algebra, dict)
  if isdefined(a, :coordinates)
    b.coordinates = Base.deepcopy_internal(a.coordinates, dict)
  end
  return b
end

################################################################################
#
#  Special elements
#
################################################################################

(O::AlgAssAbsOrd{S, T})() where {S, T} = AlgAssAbsOrdElem{S, T}(O)

one(O::AlgAssAbsOrd) = O(one(algebra(O)))

zero(O::AlgAssAbsOrd) = O(algebra(O)())

################################################################################
#
#  Element in algebra
#
################################################################################

function elem_in_algebra(x::AlgAssAbsOrdElem; copy::Bool = true) where T
  if copy
    return deepcopy(x.elem_in_algebra)
  else
    return x.elem_in_algebra
  end
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_coord(x::AlgAssAbsOrdElem)
  if x.has_coord
    return nothing
  end

  a, b = _check_elem_in_order(elem_in_algebra(x, copy = false), parent(x))
  !a && error("Not a valid order element")
  x.coordinates = b
  x.has_coord = true
  return nothing
end

################################################################################
#
#  Coordinates
#
################################################################################

function coordinates(x::AlgAssAbsOrdElem; copy::Bool = true)
  assure_has_coord(x)
  if copy
    return deepcopy(x.coordinates)
  else
    return x.coordinates
  end
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::AlgAssAbsOrdElem)
  return parent(x)(-elem_in_algebra(x, copy = false))
end

###############################################################################
#
#  Binary operations
#
###############################################################################

function *(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem)
  !check_parent(x, y) && error("Wrong parents")
  return parent(x)(elem_in_algebra(x, copy = false)*elem_in_algebra(y, copy = false))
end

function +(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem)
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)(elem_in_algebra(x, copy = false) + elem_in_algebra(y, copy = false))
  if x.has_coord && y.has_coord
    z.coordinates = [ x.coordinates[i] + y.coordinates[i] for i = 1:degree(parent(x)) ]
    z.has_coord = true
  end
  return z
end

function -(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem)
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)(elem_in_algebra(x, copy = false) - elem_in_algebra(y, copy = false))
  if x.has_coord && y.has_coord
    z.coordinates = [ x.coordinates[i] - y.coordinates[i] for i = 1:degree(parent(x)) ]
    z.has_coord = true
  end
  return z
end

function *(n::Union{Integer, fmpz}, x::AlgAssAbsOrdElem)
  O=x.parent
  y=Array{fmpz,1}(undef, O.dim)
  z=coordinates(x, copy = false)
  for i=1:O.dim
    y[i] = z[i] * n
  end
  return O(y)
end

# Computes a/b if action is :right and b\a if action is :left (and if this is possible)
function divexact(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem, action::Symbol, check::Bool = true)
  !check_parent(a, b) && error("Wrong parents")
  O = parent(a)
  c = divexact(elem_in_algebra(a, copy = false), elem_in_algebra(b, copy = false), action)
  if check
    (x, y) = _check_elem_in_order(c, O)
    !x && error("Quotient not an element of the order")
    return typeof(a)(O, c, y) # Avoid unneccessary copies
  end
  return typeof(a)(O, c)
end

divexact_right(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem, check::Bool = true) = divexact(a, b, :right, check)

divexact_left(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem, check::Bool = true) = divexact(a, b, :left, check)

################################################################################
#
#  Conversion from matrix
#
################################################################################

function elem_from_mat_row(O::AlgAssAbsOrd, M::fmpz_mat, i::Int)
  return O(fmpz[ M[i, j] for j = 1:degree(O) ])
end

function elem_to_mat_row!(M::fmpz_mat, i::Int, a::AlgAssAbsOrdElem)
  for c = 1:ncols(M)
    M[i, c] = deepcopy(coordinates(a; copy = false))[c]
  end
  return nothing
end

################################################################################
#
#  Exponentiation
#
################################################################################

function ^(x::AlgAssAbsOrdElem, y::Union{fmpz, Int})
  z = parent(x)()
  z.elem_in_algebra = elem_in_algebra(x, copy = false)^y
  return z
end

################################################################################
#
#  Equality
#
################################################################################

function ==(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem)
  if parent(a) !== parent(b)
    return false
  end
  return elem_in_algebra(a, copy = false) == elem_in_algebra(b, copy = false)
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function mul!(z::AlgAssAbsOrdElem, x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem)
  z.elem_in_algebra = mul!(elem_in_algebra(z, copy = false), elem_in_algebra(x, copy = false), elem_in_algebra(y, copy = false))
  z.has_coord = false
  return z
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdElem)
  print(io, a.elem_in_algebra)
end

################################################################################
#
#  Representation matrices
#
################################################################################

function representation_matrix(x::AlgAssAbsOrdElem, action::Symbol = :left)

  O = parent(x)
  M = basis_mat(O, copy = false)
  M1 = basis_mat_inv(O, copy = false)

  B = FakeFmpqMat(representation_matrix(elem_in_algebra(x, copy = false), action))
  B = mul!(B, M, B)
  B = mul!(B, B, M1)

  @assert B.den == 1
  return B.num
end

################################################################################
#
#  Trace
#
################################################################################

function tr(x::AlgAssAbsOrdElem)
  return FlintZZ(tr(x.elem_in_algebra))
end

function trred(x::AlgAssAbsOrdElem)
  return FlintZZ(trred(x.elem_in_algebra))
end

################################################################################
#
#  Modular exponentiation and division
#
################################################################################

function powermod(a::AlgAssAbsOrdElem, i::Union{fmpz, Int}, m::AlgAssAbsOrdIdl)
  if i < 0
    b, a = isdivisible_mod_ideal(one(parent(a)), a, m)
    @assert b "Element is not invertible modulo the ideal"
    return powermod(a, -i, m)
  end

  if i == 0
    return one(parent(a))
  end

  if i == 1
    b = mod(a, m)
    return b
  end

  if mod(i, 2) == 0
    j = div(i, 2)
    b = powermod(a, j, m)
    b = b^2
    b = mod(b, m)
    return b
  end

  b = mod(a*powermod(a, i - 1, m), m)
  return b
end

# This is mostly isdivisible in NfOrd/ResidueRing.jl
function isdivisible_mod_ideal(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl)

  iszero(y) && error("Dividing by zero")

  if iszero(x)
    return true, zero(parent(x))
  end

  O = parent(x)
  d = degree(O)
  V = zero_matrix(FlintZZ, 2*d + 1, 2*d + 1)
  V[1, 1] = fmpz(1)

  for i = 1:d
    V[1, 1 + i] = coordinates(x, copy = false)[i]
  end

  A = representation_matrix(y)
  B = basis_mat(a, copy = false)

  _copy_matrix_into_matrix(V, 2, 2, A)
  _copy_matrix_into_matrix(V, 2 + d, 2, B)

  for i = 1:d
    V[1 + i, d + 1 + i] = 1
  end

  V = hnf(V)

  for i = 2:(d + 1)
    if !iszero(V[1, i])
      return false, O()
    end
  end

  z = -O([ V[1, i] for i = (d + 2):(2*d + 1) ])
  return true, z
end

################################################################################
#
#  isone/iszero
#
################################################################################

iszero(a::AlgAssAbsOrdElem) = iszero(elem_in_algebra(a, copy = false))

isone(a::AlgAssAbsOrdElem) = isone(elem_in_algebra(a, copy = false))


text:
"""
    SEALContext

Heavyweight class that does validates encryption parameters of type `EncryptionParameters` and
pre-computes and stores several costly pre-computations.

See also: [`EncryptionParameters`](@ref)
"""
mutable struct SEALContext <: SEALObject
  handle::Ptr{Cvoid}

  function SEALContext(enc_param::EncryptionParameters;
                       expand_mod_chain=true, sec_level=SecLevelType.tc128)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:SEALContext_Create, libsealc), Clong,
                   (Ptr{Cvoid}, UInt8, Cint, Ref{Ptr{Cvoid}}),
                   enc_param, expand_mod_chain, sec_level, handleref)
    @check_return_value retval
    return SEALContext(handleref[])
  end

  function SEALContext(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::SEALContext)
  if isallocated(object)
    @check_return_value ccall((:SEALContext_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function first_parms_id(context::SEALContext)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:SEALContext_FirstParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 context, parms_id)
  @check_return_value retval
  return parms_id
end

function last_parms_id(context::SEALContext)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:SEALContext_LastParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 context, parms_id)
  @check_return_value retval
  return parms_id
end

function get_context_data(context::SEALContext, parms_id::DenseVector{UInt64})
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:SEALContext_GetContextData, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}, Ref{Ptr{Cvoid}}),
                 context, parms_id, handleref)
  @check_return_value retval
  return ContextData(handleref[], destroy_on_gc=false)
end

function key_context_data(context::SEALContext)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:SEALContext_KeyContextData, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 context, handleref)
  @check_return_value retval
  return ContextData(handleref[], destroy_on_gc=false)
end

function first_context_data(context::SEALContext)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:SEALContext_FirstContextData, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 context, handleref)
  @check_return_value retval
  return ContextData(handleref[], destroy_on_gc=false)
end

function parameter_error_message(context::SEALContext)
  len = Ref{UInt64}(0)

  # First call to obtain length (message pointer is null)
  retval = ccall((:SEALContext_ParameterErrorMessage, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt64}),
                 context, C_NULL, len)
  @check_return_value retval

  # Second call to obtain message
  message = Vector{UInt8}(undef, len[])
  retval = ccall((:SEALContext_ParameterErrorMessage, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt64}),
                 context, message, len)
  @check_return_value retval

  return String(message)
end

function using_keyswitching(context::SEALContext)
  valueref = Ref{UInt8}(0)
  retval = ccall((:SEALContext_UsingKeyswitching, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}),
                 context, valueref)
  @check_return_value retval
  return Bool(valueref[])
end

mutable struct ContextData <: SEALObject
  handle::Ptr{Cvoid}

  function ContextData(handle::Ptr{Cvoid}; destroy_on_gc=true)
    object = new(handle)
    destroy_on_gc && finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::ContextData)
  if isallocated(object)
    @check_return_value ccall((:ContextData_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function chain_index(context_data::ContextData)
  index = Ref{UInt64}(0)
  retval = ccall((:ContextData_ChainIndex, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 context_data, index)
  @check_return_value retval
  return Int(index[])
end

function parms(context_data::ContextData)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:ContextData_Parms, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 context_data, handleref)
  @check_return_value retval
  return EncryptionParameters(handleref[])
end

function parms_id(context_data::ContextData)
  enc_parms = parms(context_data)
  return parms_id(enc_parms)
end

function total_coeff_modulus_bit_count(context_data::ContextData)
  bit_count = Ref{Cint}(0)
  retval = ccall((:ContextData_TotalCoeffModulusBitCount, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cint}),
                 context_data, bit_count)
  @check_return_value retval
  return Int(bit_count[])
end

function qualifiers(context_data::ContextData)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:ContextData_Qualifiers, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 context_data, handleref)
  @check_return_value retval
  return EncryptionParameterQualifiers(handleref[], destroy_on_gc=false)
end

function next_context_data(context_data::ContextData)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:ContextData_NextContextData, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 context_data, handleref)
  @check_return_value retval
  if handleref[] == C_NULL
    return nothing
  else
    return ContextData(handleref[], destroy_on_gc=false)
  end
end

mutable struct EncryptionParameterQualifiers <: SEALObject
  handle::Ptr{Cvoid}

  function EncryptionParameterQualifiers(handle::Ptr{Cvoid}; destroy_on_gc=true)
    object = new(handle)
    destroy_on_gc && finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::EncryptionParameterQualifiers)
  if isallocated(object)
    @check_return_value ccall((:EncryptionParameterQualifiers_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function using_batching(epq::EncryptionParameterQualifiers)
  valueref = Ref{UInt8}(0)
  retval = ccall((:EPQ_UsingBatching, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}),
                 epq, valueref)
  @check_return_value retval
  return Bool(valueref[])
end


text:
using Documenter

push!(LOAD_PATH, "../src/")
using XPA

DEPLOYDOCS = (get(ENV, "CI", nothing) == "true")

makedocs(
    sitename = "XPA.jl Package",
    format = Documenter.HTML(
        prettyurls = DEPLOYDOCS,
    ),
    authors = "Éric Thiébaut and contributors",
    pages = ["index.md", "install.md", "intro.md",
             "client.md", "server.md", "misc.md", "library.md"]
)

if DEPLOYDOCS
    deploydocs(
        repo = "github.com/emmt/XPA.jl.git",
    )
end


text:
using LiveServer, Test, Crayons, Sockets, HTTP
const LS = LiveServer
# depending on the OS, the resolution of stat(f).mtime can be over a second.
const FS_WAIT = 1.1

include("utils.jl")
include("file_watching.jl")
include("server.jl")


text:
@cache mutable struct ISSEMCache{uType,rateType,JType,WType,JC,UF,
                          N,noiseRateType,F,randType} <:
                          StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  tmp::uType
  gtmp::noiseRateType
  gtmp2::rateType
  J::JType
  W::WType
  jac_config::JC
  linsolve::F
  uf::UF
  nlsolve::N
  dW_cache::randType
end

function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  @iipnlcachefields
  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    gtmp2 = gtmp
    dW_cache = nothing
  else
    gtmp2 = zero(rate_prototype)
    dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  @iipnlsolve
  ISSEMCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,J,W,jac_config,linsolve,uf,
                  nlsolve,dW_cache)
end

mutable struct ISSEMConstantCache{F,N} <: StochasticDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
  @oopnlcachefields
  γ, c = alg.theta,zero(t)
  @oopnlsolve
  ISSEMConstantCache(uf,nlsolve)
end

@cache mutable struct ISSEulerHeunCache{uType,rateType,JType,WType,JC,UF,N,
                                 noiseRateType,F,randType} <:
                                 StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  tmp::uType
  gtmp::noiseRateType
  gtmp2::rateType
  gtmp3::noiseRateType
  J::JType
  W::WType
  jac_config::JC
  linsolve::F
  uf::UF
  nlsolve::N
  dW_cache::randType
end

function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  @iipnlcachefields

  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)

  if is_diagonal_noise(prob)
      gtmp3 = gtmp2
      dW_cache = nothing
  else
      gtmp3 = zero(noise_rate_prototype)
      dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  @iipnlsolve
  ISSEulerHeunCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,gtmp3,
                         J,W,jac_config,linsolve,uf,nlsolve,dW_cache)
end

mutable struct ISSEulerHeunConstantCache{F,N} <: StochasticDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
  @oopnlcachefields
  γ, c = alg.theta,zero(t)
  @oopnlsolve
  ISSEulerHeunConstantCache(uf,nlsolve)
end


text:
using Descartes
using GeometryBasics: Mesh

function beam(beam_size = [50,10,10],
              hole_ct = 5,
              hole_d = 3)

    # computations
    hole_interval = beam_size[1]/(hole_ct + 1)

    c = Square(beam_size)

    for i = 1:hole_ct
        h = translate([hole_interval*i, beam_size[2]/2])Circle(hole_d/2)
        c = diff(c, h)
    end
    LinearExtrude(c, beam_size[3])
end

m = Mesh(beam())

#save("2d_beam.ply",m)


text:
module CUDA

using CuArrays

CuArrays.cudnn_available() && include("cudnn.jl")

end


text:
using Test
using LinearAlgebra
using IntervalArithmetic
using PolynomialBasis
# using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function allapprox(v1, v2; tol = 1e-14)
    @assert length(v1) == length(v2)
    flags = [isapprox(v1[i], v2[i], atol = tol) for i = 1:length(v1)]
    return all(flags)
end

function circle_distance_function(coords, center, radius)
    npts = size(coords)[2]
    return [norm(coords[:, i] - center) - radius for i = 1:npts]
end

function plane_distance_function(coords, normal, x0)
    return (coords .- x0)' * normal
end

function plot_zero_levelset(poly)
    x = -1:1e-2:1
    contour(x, x, (x, y) -> poly(x, y), levels = [0.0])
end

struct FunctionAndGradient
    func::Any
    grad::Any
    perturbation::Any
end

function (f::FunctionAndGradient)(x)
    return f.func(x) + perturbation
end

function PolynomialBasis.gradient(f::FunctionAndGradient, x)
    return f.grad(x)
end

#######################################################################
# Test a curved interface and subdivision
radius = 0.5
center = [0.0, 0.0]
poly = InterpolatingPolynomial(1, 2, 3)
dpoly = InterpolatingPolynomial(2, 2, 3)
coeffs = circle_distance_function(poly.basis.points, center, radius)
update!(poly, coeffs)
update_interpolating_gradient!(dpoly,poly)
testarea = 4.0 - pi * 0.5^2

xL,xR = [-1.0,-1.0],[1.,1.]

numqp = 5
quad = IDQ.area_quadrature(poly,dpoly, +1, xL,xR, numqp)
@test isapprox(sum(quad.weights), testarea, atol = 1e-1)

squad = IDQ.surface_quadrature(poly,dpoly, xL,xR, numqp)
@test isapprox(sum(quad.weights), pi, atol = 1e-1)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
using Distributed
using ProgressMeter
using Logging
using SharedArrays
using JLD2
using Dates
using Parallelism
# using Config

include("slurm.jl")
using .ClusterManagers

IN_SLURM() = ("SLURM_JOBID" ∈ keys(ENV)) && ("SLURM_NTASKS" ∈ keys(ENV))

"""
    job(experiment::Experiment; kwargs...)
    job(experiment_file, exp_dir, args_iter; kwargs...)
    job(experiment::Experiment, job_id; kwargs...)
    job(experiment_file, exp_dir, args_iter, job_id; kwargs...)

Run a job specified by the experiment.
"""
job(exp::Experiment; kwargs...) =
    job(exp.file,
        exp.dir,
        exp.args_iter;
        exp_module_name=exp.module_name,
        exp_func_name=exp.func_name,
        exception_dir="$(exp.dir)/except/exp_0x$(string(exp.hash, base=16))",
        checkpoint_name="$(exp.dir)/checkpoints/exp_0x$(string(exp.hash, base=16))",
        kwargs...)

function job(experiment_file,
             exp_dir,
             args_iter;
             kwargs...)
    if "SLURM_ARRAY_TASK_ID" in keys(ENV)
        @info "This is an array Job! Time to get task and start job."
        task_id = parse(Int64, ENV["SLURM_ARRAY_TASK_ID"])
        @time task_job(experiment_file,
                       exp_dir,
                       args_iter,
                       task_id;
                       kwargs...)
    else
        @time parallel_job(experiment_file,
                           exp_dir,
                           args_iter;
                           kwargs...)
    end
end

function job(experiment_file::AbstractString,
             exp_dir::AbstractString,
             args_iter,
             job_id;
             kwargs...)
    @info "This is a task job! ID is $(job_id)"
    @time task_job(experiment_file,
                   exp_dir,
                   args_iter,
                   job_id;
                   kwargs...)
end

function create_procs(num_workers, project, job_file_dir)
    # assume started fresh julia instance...

    pids = Array{Int64, 1}()
    exc_opts = Base.JLOptions()
    color_opt = "no"
    if exc_opts.color == 1
        color_opt = "yes"
    end

    if IN_SLURM()
        if nworkers() == 1
            num_workers = parse(Int64, ENV["SLURM_NTASKS"])
            pids = addprocs(SlurmManager(num_workers);
                            exeflags=["--project=$(project)", "--color=$(color_opt)"],
                            job_file_loc=job_file_dir)
        else
            pids = procs()
        end
    else
        if nworkers() == 1
            pids = addprocs(num_workers;
                            exeflags=["--project=$(project)", "--color=$(color_opt)"])
        elseif nworkers() < num_workers
            pids = addprocs((num_workers) - nworkers();
                            exeflags=["--project=$(project)", "--color=$(color_opt)"])
        else
            pids = procs()
        end
    end
    return fetch(pids)
end

function run_experiment(exp_func, job_id, args, extra_args, exception_loc;
                         expand_args=false,
                         verbose=false,
                         store_exceptions=true,
                         skip_exceptions=false)

    run_exp = !(isfile(joinpath(exception_loc, "job_$(job_id).exc")))

    if run_exp || !skip_exceptions
        try
            if expand_args
                exp_func(args..., extra_args...)
            else
                exp_func(args, extra_args...)
            end
            return true
        catch ex
            if isa(ex, InterruptException)
                throw(InterruptException())
            end
            if verbose
                @warn "Exception encountered for job: $(job_id)"
            end
            if store_exceptions
                exception_file(
                    joinpath(exception_loc, "job_$(job_id).exc"),
                    job_id, ex, stacktrace(catch_backtrace()))
            end
            return false
        end
    elseif verbose
        @warn  "Not running job for $(job_id)"
    end

    return true

end

"""
    parallel_job

Run a parallel job over the arguments presented by args_iter. `args_iter` can be a enumeration OR ArgIterator. Each job will be dedicated to a specific
task. The experiment *must* save its own data! As this is not handled by this function (although could be added in the future.)
"""
function parallel_job(experiment_file::AbstractString,
                      exp_dir::AbstractString,
                      args_iter;
                      exp_module_name::Union{String, Symbol}=:Main,
                      exp_func_name::Union{String, Symbol}=:main_experiment,
                      num_workers=Sys.CPU_THREADS - 1,
                      project=".",
                      extra_args=[],
                      exception_dir="except",
                      job_file_dir="",
                      checkpoint_name="",
                      store_exceptions=true,
                      verbose=false,
                      skip_exceptions=false,
                      expand_args=false)

    #######
    #
    # Preamble - Add processes, initialized shared memory
    #
    ######
    if job_file_dir == ""
        job_file_dir = joinpath(exp_dir, "jobs")
    end

    if exception_dir == "except"
        exception_dir = joinpath(exp_dir, exception_dir)
    end
    if store_exceptions && !isdir(exception_dir)
        mkpath(exception_dir)
    end

    checkpointing = true
    checkpoint_file = checkpoint_name
    if checkpoint_name == ""
        checkpointing = false
    elseif splitext(checkpoint_file)[end] != ".jld2"
        checkpoint_file = checkpoint_file * ".jld2"
    end

    n = length(args_iter)

    finished_jobs_arr = fill(false, n)
    if checkpointing
        if !isdir(dirname(checkpoint_file))
            mkpath(dirname(checkpoint_file))
        end
        if isfile(checkpoint_file)
            JLD2.@load checkpoint_file finished_jobs_arr
        else
            JLD2.@save checkpoint_file finished_jobs_arr
        end
    end
    done_jobs = finished_jobs_arr

    @info "Number of Jobs left: $(n - sum(done_jobs))/$(n)"

    if all(done_jobs)
        @info "All jobs finished!"
        return findall((x)->x==false, finished_jobs_arr)
    end

    if !(experiment_file[1] == '/' || experiment_file[1] == '\\')
        experiment_file = joinpath(pwd(), experiment_file)
    end

    #########
    #
    # Meaty middle: Compiling code, running jobs, managing which jobs fail.
    #
    ########

    # job_id_channel: a job id will appear here if a job is finished.
    job_id_channel = RemoteChannel(()->Channel{Int}(min(1000, length(args_iter))), 1)
    prg_channel = RemoteChannel(()->Channel{Bool}(min(1000, length(args_iter))), 1)

    # Include on first proc for pre-compiliation
    @info "pre-compile"
    @everywhere begin
        include($experiment_file)
    end

    # create processes.
    pids = create_procs(num_workers, project, job_file_dir)
    println(nworkers(), " ", pids)

    try

        mod_str = string(exp_module_name)
        func_str = string(exp_func_name)

        @everywhere const global RP_exp_file=$experiment_file

        @everywhere begin
            eval(:(using Reproduce))
            eval(:(using Distributed))
            eval(:(using SharedArrays))

            include(RP_exp_file)
            @info "$(RP_exp_file) included on process $(myid())"
            mod = $mod_str=="Main" ? Main : getfield(Main, Symbol($mod_str))
            const global RP_exp_func = getfield(mod, Symbol($func_str))
            @info "Experiment built on process $(myid())"
        end

        pgm = ProgressMeter.Progress(length(args_iter))

        @sync begin
            @async while Distributed.take!(prg_channel)
                ProgressMeter.next!(pgm)
                JLD2.@load checkpoint_file finished_jobs_arr
                while isready(job_id_channel)
                    new_job_id = take!(job_id_channel)
                    finished_jobs_arr[new_job_id] = true
                end
                JLD2.@save checkpoint_file finished_jobs_arr
            end

            @sync begin
                robust_pmap(args_iter) do (job_id, args)
                    if !checkpointing || !done_jobs[job_id]
                        finished = run_experiment(Main.RP_exp_func, job_id, args, extra_args, exception_dir;
                                                  expand_args=expand_args,
                                                  verbose=verbose,
                                                  store_exceptions=store_exceptions,
                                                  skip_exceptions=skip_exceptions)
                        if finished
                            Distributed.put!(job_id_channel, job_id)
                        end
                    end
                    Distributed.put!(prg_channel, true)
                    yield()
                end
                Distributed.put!(prg_channel, false)
            end
        end

        if checkpointing
            JLD2.@load checkpoint_file finished_jobs_arr
            return findall((x)->x==false, finished_jobs_arr)
        else
            finished_jobs_bool = fill(false, n)
            while isready(job_id_channel)
                job_id = take!(job_id_channel)
                finished_jobs_bool[job_id] = true
            end
            return findall((x)->x==false, finished_jobs_bool)
        end
    catch ex
        Distributed.interrupt()
        rethrow()
    end
end

function task_job(experiment_file::AbstractString,
                  exp_dir::AbstractString,
                  args_iter, task_id::Integer;
                  exp_module_name::Union{String, Symbol}=:Main,
                  exp_func_name::Union{String, Symbol}=:main_experiment,
                  project=".",
                  extra_args=[],
                  exception_dir="except",
                  checkpoint_name="",
                  store_exceptions=true,
                  verbose=false,
                  skip_exceptions=false,
                  expand_args=false)

    if exception_dir == "except"
        exception_dir = joinpath(exp_dir, exception_dir)
    end
    if store_exceptions && !isdir(exception_dir)
        _safe_mkpath(exception_dir)
    end

    checkpointing = false
    checkpoint_folder = checkpoint_name
    if checkpoint_folder != ""
        _safe_mkpath(checkpoint_folder)
        checkpointing = true
        open(joinpath(checkpoint_folder, "job_$(task_id)"), "w") do f
            write(f, "Not Done")
        end
    end

    mod_str = string(exp_module_name)
    func_str = string(exp_func_name)

    @everywhere begin
        include($experiment_file)
        mod = $mod_str=="Main" ? Main : getfield(Main, Symbol($mod_str))
        const global exp_func = getfield(mod, Symbol($func_str))
    end

    job_id = task_id
    args = collect(args_iter)[task_id][2]
    ret = @sync @async begin
        run_experiment(Main.exp_func, job_id, args, extra_args, exception_dir;
                       expand_args=expand_args,
                       verbose=verbose,
                       store_exceptions=store_exceptions,
                       skip_exceptions=skip_exceptions)
    end

    if checkpointing
        rm(joinpath(checkpoint_folder, "job_$(task_id)"))
        if ret == false
            open(joinpath(checkpoint_folder, "job_$(task_id)"), "w") do f
                write(f, "exception while running job! Check exceptions for more details.")
            end
        end
    end

    return true
end