text:
# ==============================================================================
# Check type of Max-Plus and type promotions on scalars
# ==============================================================================

# ==============================================================================
# Scalar constructor
# ==============================================================================

# Min-Plus

a = MI(1.0)
@test typeof(a) == MI
@test a isa MI
@test a.v isa Float64
@test a.v == 1.0
@test isnan(a) == false

b = MI(1)
@test typeof(b) == MI
@test b isa MI
@test b.v isa Float64
@test b.v == 1
@test isnan(b) == false

c = MI(-1)
@test c isa MI
@test c.v isa Float64
@test c.v == -1.0
@test isnan(c) == false

d = -MI(1)
@test d isa MI
@test d.v isa Float64
@test d.v == -1.0
@test isnan(d) == false

e = MI(NaN)
@test e isa MI
@test e.v isa Float64
@test isnan(e) == true

f = Trop(MI, NaN)
@test f isa MI
@test f.v isa Float64
@test isnan(f) == false
@test f.v == Inf

g = Trop(MI, 42.0)
@test g isa MI
@test g.v isa Float64
@test isnan(g) == false
@test g.v == 42.0

# Max-Plus

a = MP(1.0)
@test typeof(a) == MP
@test a isa MP
@test a.v isa Float64
@test a.v == 1.0
@test isnan(a) == false

b = MP(1)
@test typeof(b) == MP
@test b isa MP
@test b.v isa Float64
@test b.v == 1
@test isnan(b) == false

c = MP(-1)
@test c isa MP
@test c.v isa Float64
@test c.v == -1.0
@test isnan(c) == false

d = -MP(1)
@test d isa MP
@test d.v isa Float64
@test d.v == -1.0
@test isnan(d) == false

e = MP(NaN)
@test e isa MP
@test e.v isa Float64
@test isnan(e) == true

f = Trop(MP, NaN)
@test f isa MP
@test f.v isa Float64
@test isnan(f) == false
@test f.v == -Inf

g = Trop(MP, 42.0)
@test g isa MP
@test g.v isa Float64
@test isnan(g) == false
@test g.v == 42.0

# ==============================================================================
# Scalar copy constructor
# ==============================================================================

a = MP(MP(1))
@test a isa MP
@test a.v == 1.0

b = MI(MI(2))
@test b isa MI
@test b.v == 2.0

c = MP(MI(3))
@test c isa MP
@test c.v == 3.0

d = MI(MP(4))
@test d isa MI
@test d.v == 4.0

# ==============================================================================
# Boolean constructor
# ==============================================================================

@test MP(true) == mp1
@test MP(false) == mp0
@test MI(true) == mi1
@test MI(false) == mi0

# ==============================================================================
# Signs
# ==============================================================================

@test sign(MP(0)) == 0
@test sign(MP(5)) == 1
@test sign(MP(-5)) == -1
@test sign(MI(0)) == 0
@test sign(MI(5)) == 1
@test sign(MI(-5)) == -1
@test sign(mp1) == 0
@test sign(mp0) == -1
@test sign(mi1) == 0
@test sign(mi0) == 1

# ==============================================================================
# Scalar comparaisons
# ==============================================================================

# Min-Plus scalars

a = MI(3.0)
@test (a == a) == true
@test (a != a) == (a ≠ a) == false
@test (a >= a) == (a ≥ a) == true
@test (a <= a) == (a ≤ a) == true
@test (a > a) == false
@test (a < a) == false

# Min-Plus vs. non Max-Plus comparaison

@test (a == 3.0) == true
@test (a == 4.0) == false
@test (a != 4.0) == true
@test (a < 4.0) == true
@test (a <= 4.0) == true
@test (a > 4.0) == false
@test (a >= 4.0) == false
@test (4.0 > a) == true
@test (4.0 >= a) == true

# Max-Plus scalars

a = MP(3.0)
@test (a == a) == true
@test (a != a) == (a ≠ a) == false
@test (a >= a) == (a ≥ a) == true
@test (a <= a) == (a ≤ a) == true
@test (a > a) == false
@test (a < a) == false

# Max-Plus vs. non Max-Plus comparaison

@test (a == 3.0) == true
@test (a == 4.0) == false
@test (a != 4.0) == true
@test (a < 4.0) == true
@test (a <= 4.0) == true
@test (a > 4.0) == false
@test (a >= 4.0) == false
@test (4.0 > a) == true
@test (4.0 >= a) == true

# ==============================================================================
# Max/Min-Plus scalar to classic algebra scalar conversion
# ==============================================================================

# Min-Plus

@test plustimes(MI(3.0)) isa Float64
@test plustimes(MI(3.0)) == 3.0

@test plustimes(MI(-3.0)) isa Float64
@test plustimes(MI(-3.0)) == -3.0

# Max-Plus

@test plustimes(MP(2.0)) isa Float64
@test plustimes(MP(2.0)) == 2.0

@test plustimes(MP(-2.0)) isa Float64
@test plustimes(MP(-2.0)) == -2.0

# ==============================================================================
# Neutral / absorbing constant elements
# ==============================================================================

# Min-Plus

@test mi0 isa MI
@test mi0.v isa Float64
@test mi0.v == Inf
@test iszero(mp0) == true
@test isone(mp0) == false

@test mi1 isa MI
@test mi1.v isa Float64
@test mi1.v == 0.0
@test iszero(mi1) == false
@test isone(mi1) == true

@test mitop isa MI
@test mitop.v isa Float64
@test mitop.v == -Inf
@test iszero(mitop) == false
@test isone(mitop) == false

# Max-Plus

@test mp0 isa MP
@test mp0.v isa Float64
@test mp0.v == -Inf
@test iszero(mp0) == true
@test isone(mp0) == false

@test mp1 isa MP
@test mp1.v isa Float64
@test mp1.v == 0.0
@test iszero(mp1) == false
@test isone(mp1) == true

@test mptop isa MP
@test mptop.v isa Float64
@test mptop.v == Inf
@test iszero(mptop) == false
@test isone(mptop) == false

# ==============================================================================
# zero / one elements
# ==============================================================================

# Min-Plus

@test zero(MI) == MI(Inf) == Inf
@test zero(MI) == mi0 == Inf
@test zero(MI(42.0)) == mi0 == Inf
@test zero(MI(42)) == mi0 == Inf

@test one(MI) == MI(0.0) == 0.0
@test one(MI) == mi1 == 0.0
@test one(MI(42.0)) == mi1 == 0.0
@test one(MI) == mi1 == 0.0
@test one(MI(42)) == mi1 == 0

# Max-Plus

@test zero(MP) == MP(-Inf) == -Inf
@test zero(MP) == mp0 == -Inf
@test zero(MP(42.0)) == mp0 == -Inf
@test zero(MP(42)) == mp0 == -Inf

@test one(MP) == MP(0.0) == 0.0
@test one(MP) == mp1 == 0.0
@test one(MP(42.0)) == mp1 == 0.0
@test one(MP) == mp1 == 0.0
@test one(MP(42)) == mp1 == 0

# ==============================================================================
# Operations on scalars and neutral elements and commutativity
# ==============================================================================

# Min-Plus

b = MI(3.0)
@test b + mi0 == mi0 + b == MI(min(3.0, Inf)) == MI(min(Inf, 3.0)) == b == MI(3.0)
@test b * mi1 == mi1 * b == MI(3.0 + 0.0) == MI(0.0 + 3.0) == b == MI(3.0)

@test b * b == MI(3.0 + 3.0) == MI(6.0)
@test b + b == MI(min(3.0, 3.0)) == MI(3.0)

a = MI(5.0)
b = MI(3.0)
@test a + b == MI(min(5.0, 3.0)) == b + a == MI(min(3.0, 5.0)) == MI(3.0)
@test a * b == MI(5.0 + 3.0) == b * a == MI(3.0 + 5.0) == MI(8.0)

# Max-Plus

b = MP(3.0)
@test b + mp0 == mp0 + b == MP(max(3.0, -Inf)) == MP(max(-Inf, 3.0)) == b == MP(3.0)
@test b * mp1 == mp1 * b == MP(3.0 + 0.0) == MP(0.0 + 3.0) == b == MP(3.0)

@test b * b == MP(3.0 + 3.0) == MP(6.0)
@test b + b == MP(max(3.0, 3.0)) == MP(3.0)

a = MP(5.0)
b = MP(3.0)
@test a + b == MP(max(5.0, 3.0)) == b + a == MP(max(3.0, 5.0)) == MP(5.0)
@test a * b == MP(5.0 + 3.0) == b * a == MP(3.0 + 5.0) == MP(8.0)

# ==============================================================================
# Operations on scalars of different types
# ==============================================================================

# Min-Plus

b = MI(3.0)
@test b + 3.0 == 3.0 + b == b == MI(min(3.0, 3.0)) == MI(3.0)
@test b + 3  == 3 + b == b == MI(min(3.0, 3.0)) == MI(3.0)
@test b * 3.0 == 3.0 * b == MI(3.0 + 3.0) == MI(6.0)
@test b * 3  == 3 * b == MI(3.0 + 3.0) == MI(6.0)

# Max-Plus

b = MP(3.0)
@test b + 3.0 == 3.0 + b == b == MP(max(3.0, 3.0)) == MP(3.0)
@test b + 3  == 3 + b == b == MP(max(3.0, 3.0)) == MP(3.0)
@test b * 3.0 == 3.0 * b == MP(3.0 + 3.0) == MP(6.0)
@test b * 3  == 3 * b == MP(3.0 + 3.0) == MP(6.0)

# ==============================================================================
#  Distributivity of operations
# ==============================================================================

# Min-Plus

a = MI(5.0)
b = MI(3.0)
c = MI(1.0)
@test (a + b) + c == a + (b + c) == a + b + c == MI(min(1.0, 3.0, 5.0)) == MI(1.0)
@test (a * b) * c == a * (b * c) == a * b * c == MI(1.0 + 3.0 + 5.0) == MI(9.0)
@test (a + b) * c == MI(min(5.0, 3.0) + 1.0) == MI(4.0)
@test (a * c) + (b * c) == MI(min(5.0 + 1.0, 3.0 + 1.0)) == MI(4.0)

# Max-Plus

a = MP(5.0)
b = MP(3.0)
c = MP(1.0)
@test (a + b) + c == a + (b + c) == a + b + c == MP(max(1.0, 3.0, 5.0)) == MP(5.0)
@test (a * b) * c == a * (b * c) == a * b * c == MP(1.0 + 3.0 + 5.0) == MP(9.0)
@test (a + b) * c == MP(max(5.0, 3.0) + 1.0) == MP(6.0)
@test (a * c) + (b * c) == MP(max(5.0 + 1.0, 3.0 + 1.0)) == MP(6.0)

# ==============================================================================
# Residuation
# ==============================================================================

# Min-Plus TODO

# Max-Plus

@test MP(3) / MP(3) == -(MP(3) * MP(-3)) == MP(0.0)
@test MP(3) / MP(5) == -(MP(5.0) * MP(-3.0)) == MP(-2.0)
@test MP(3) \ MP(3) == -(MP(3) * MP(-3)) == MP(0.0)
@test MP(3) \ MP(6) == -(MP(3.0) * MP(-6.0)) == MP(3.0)
@test MP(3) \ mp0 == mp0
@test MP(3) \ mp1 == MP(-3)

# ==============================================================================
# Forbiden minus operator
# ==============================================================================

@test_throws ErrorException("Minus operator does not exist in min+ algebra") MI(3) - MI(3)
@test_throws ErrorException("Minus operator does not exist in min+ algebra") MI(3) - 3
@test_throws ErrorException("Minus operator does not exist in min+ algebra") 3 - MI(3)
@test_throws ErrorException("Minus operator does not exist in max+ algebra") MP(3) - MP(3)
@test_throws ErrorException("Minus operator does not exist in max+ algebra") MP(3) - 3
@test_throws ErrorException("Minus operator does not exist in max+ algebra") 3 - MP(3)

# ==============================================================================
# Scalar power operator
# ==============================================================================

# Min-Plus TODO

# Max-Plus

@test MP(2)^4  == MP(2 * 4) == MP(8)
@test MP(0)^0  == MP(0 * 0) == MP(0)
@test MP(2)^0  == MP(2 * 0) == MP(0)
@test MP(2)^-3 == MP(2)^(-3) == MP(2 * -3) == MP(-6)
@test MP(2)^0.5 == MP(2 * 0.5) == MP(1.0)
@test MP(2)^(-0.5) == MP(2 * -0.5) == MP(-1.0)
@test mp0^0 == MP(0.0)
@test mp0^2 == mp0
@test_broken mp0^(-2) == mp0
@test inv(MP(5)) == MP(5)^-1 == MP(-5)

# ==============================================================================
# Max-Plus min operator
# ==============================================================================

@test min(MP(3.0), mp0) == min(mp0, MP(3.0)) == mp0
@test min(MP(3.0), mp1) == min(mp1, MP(3.0)) == mp1
@test min(MP(1), MP(2)) == min(1, MP(2)) == min(MP(1), 2) == MP(1)
#@test min(MP([10 1; 10 1]), MP([4 5; 6 5])) == MP([4 1; 6 1])
#@test min(sparse(MP([10.0 1; 10 1])), sparse(MP([4.0 5; 6 5]))) == sparse(MP([4.0 1; 6 1]))
#@test min(sparse(MP([10 1; 10 1])), sparse(MP([4 5; 6 5]))) == sparse(MP([4 1; 6 1]))
#@test min(sparse(eye(MP,2,2)), spzeros(MP,2,2)) == spzeros(MP,2,2)
# @test min(eye(MP,2,2), ones(2,2)) == eye(MP,2,2)

# ==============================================================================
# Border cases operations
# ==============================================================================

# Min-Plus

@test mi0 + mi0 == mi0
@test mi0 + mi1 == 0
@test mi0 + mitop == mitop
@test mi1 + mi0 == mi1
@test mi1 + mi1 == mi1
@test mi1 + mitop == mitop
@test mitop + mi0 == mitop
@test mitop + mi1 == mitop
@test mitop + mitop == mitop

@test mi0 * mi0 == mi0
@test mi0 * mi1 == mi0
#@test mi0 * mitop == mitop
@test mi1 * mi0 == mi0
@test mi1 * mi1 == mi1
@test mi1 * mitop == mitop
#@test mitop * mi0 == mitop
@test mitop * mi1 == mitop
@test mitop * mitop == mitop

#@test mi0 / mi0 == mi0
#@test mi0 / mi1 == mi0
#@test mi0 / mitop == mi0
#@test mi1 / mi0 == mi1
#@test mi1 / mi1 == mi1
#@test mi1 / mitop == mi1
#@test mitop / mi0 == mi1
#@test mitop / mi1 == mi1
#@test mitop / mitop == mi1

#@test mi0 \ mi0 == mi0
#@test mi0 \ mi1 == mi0
#@test mi0 \ mitop == mi0
#@test mi1 \ mi0 == mi1
#@test mi1 \ mi1 == mi1
#@test mi1 \ mitop == mi1
#@test mitop \ mi0 == mi1
#@test mitop \ mi1 == mi1
#@test mitop \ mitop == mi1

# Max-Plus

@test mp0 + mp0 == mp0
@test mp0 + mp1 == 0
@test mp0 + mptop == mptop
@test mp1 + mp0 == mp1
@test mp1 + mp1 == mp1
@test mp1 + mptop == mptop
@test mptop + mp0 == mptop
@test mptop + mp1 == mptop
@test mptop + mptop == mptop

@test mp0 * mp0 == mp0
@test mp0 * mp1 == mp0
@test mp0 * mptop == mp0
@test mp1 * mp0 == mp0
@test mp1 * mp1 == mp1
@test mp1 * mptop == mptop
@test mptop * mp0 == mp0
@test mptop * mp1 == mptop
@test mptop * mptop == mptop

@test mp0 / mp0 == mptop
@test mp0 / mp1 == mp0
@test mp0 / mptop == mp0
@test mp1 / mp0 == mptop
@test mp1 / mp1 == mp1
@test mp1 / mptop == mp0
@test mptop / mp0 == mptop
@test mptop / mp1 == mptop
@test mptop / mptop == mptop

@test mp0 \ mp0 == mptop
@test mp0 \ mp1 == mptop
@test mp0 \ mptop == mptop
@test mp1 \ mp0 == mp0
@test mp1 \ mp1 == mp1
@test mp1 \ mptop == mptop
@test mptop \ mp0 == mp0
@test mptop \ mp1 == mp0
@test mptop \ mptop == mptop

# ==============================================================================
# Other operations
# ==============================================================================

# Min-Plus

@test abs2(MI(3.0)) == MI(6.0)
@test abs(MI(-3.0)) == MI(3.0)
@test abs(mi0) == mi0
@test float(MI(2)) == MI(2)
#@test log(MI(2)) == MI(log(2.0)) == log(2.0)
#@test round(MI(1.7)) == MI(round(1.7)) == round(1.7)
#@test floor(MI(1.7)) == MI(floor(1.7)) == floor(1.7)

# Max-Plus

@test abs2(MP(3.0)) == MP(6.0)
@test abs(MP(-3.0)) == MP(3.0)
@test abs(mp0) == mptop
@test float(MP(2)) == MP(2)
#@test log(MP(2)) == MP(log(2.0)) == log(2.0)
#@test round(MP(1.7)) == MP(round(1.7)) == round(1.7)
#@test floor(MP(1.7)) == MP(floor(1.7)) == floor(1.7)

# ==============================================================================
# Max-Plus matrix comparaisons
# ==============================================================================

# Dense Max-Plus matrix comparaison

B = [MP(1) MP(2); MP(3) MP(4)]
@test (B == B) == true
@test (B != B) == (B ≠ B) == false
@test (B .>= B) == (B .≥ B) == [true true; true true]
@test (B .<= B) == (B .≤ B) == [true true; true true]
@test (B .> B) == [false false; false false]
@test (B .< B) == [false false; false false]

# TODO Dense vector comparaison

# Sparse Max-Plus matrix comparaison

S = sparse([MP(1) MP(2); MP(3) MP(4)])
@test (B == B) == true
@test (B != B) == (B ≠ B) == false
@test (B .>= B) == (B .≥ B) == [true true; true true]
@test (B .<= B) == (B .≤ B) == [true true; true true]
@test (B .> B) == [false false; false false]
@test (B .< B) == [false false; false false]

# TODO Sparse vector comparaison

# Sparse/Dense Max-Plus matrix comparaison

@test (S == B) == (B == S) == true

# Sort is using isless()
v = MP([3, 1, 2]);
@test sort(v) == MP([1, 2, 3])

# ==============================================================================
# Min-Plus matrix comparaisons
# ==============================================================================

# Dense Min-Plus matrix comparaison

B = [MI(1) MI(2); MI(3) MI(4)]
@test (B == B) == true
@test (B != B) == (B ≠ B) == false
@test (B .>= B) == (B .≥ B) == [true true; true true]
@test (B .<= B) == (B .≤ B) == [true true; true true]
@test (B .> B) == [false false; false false]
@test (B .< B) == [false false; false false]

# Sparse Min-Plus matrix comparaison

S = sparse([MI(1) MI(2); MI(3) MI(4)])
@test (B == B) == true
@test (B != B) == (B ≠ B) == false
@test (B .>= B) == (B .≥ B) == [true true; true true]
@test (B .<= B) == (B .≤ B) == [true true; true true]
@test (B .> B) == [false false; false false]
@test (B .< B) == [false false; false false]

# Sparse/Dense Min-Plus matrix comparaison

@test (S == B) == (B == S) == true

# Sort is using isless()
v = MI([3, 1, 2]);
@test sort(v) == MI([1, 2, 3])

# ==============================================================================
# Range construction
# ==============================================================================

# Min-Plus

A = MI(1:3)
@test A == [MI(1); MI(2); MI(3)]

B = MI(1.0:0.5:3.0)
@test B == [MI(1.0); MI(1.5); MI(2.0); MI(2.5); MI(3.0)]

# Max-Plus

A = MP(1:3)
@test A == [MP(1); MP(2); MP(3)]

B = MP(1.0:0.5:3.0)
@test B == [MP(1.0); MP(1.5); MP(2.0); MP(2.5); MP(3.0)]

# ==============================================================================
# Dense matrix constructor and type promotion/contamination
# ==============================================================================

# Min-Plus

A = [MI(1) MI(2); MI(3) MI(4)]
@test A isa Matrix{MI}

B = MI([1 2; 3 4])
@test B isa Matrix{MI}

C = [MI(1) 2; 3 4]
@test C isa Matrix{MI}

D = [MI(1.0) 2; 3 4]
@test D isa Matrix{MI}

E = [MI(1) 2.0; 3 4]
@test E isa Matrix{MI}

F = [1 MI(2.0); 3 4]
@test F isa Matrix{MI}

@test A == B == C == D == E == F

# Max-Plus

A = [MP(1) MP(2); MP(3) MP(4)]
@test A isa Matrix{MP}

B = MP([1 2; 3 4])
@test B isa Matrix{MP}

C = [MP(1) 2; 3 4]
@test C isa Matrix{MP}

D = [MP(1.0) 2; 3 4]
@test D isa Matrix{MP}

E = [MP(1) 2.0; 3 4]
@test E isa Matrix{MP}

F = [1 MP(2.0); 3 4]
@test F isa Matrix{MP}

@test A == B == C == D == E == F

# ==============================================================================
# Max-Plus sparse constructor
# ==============================================================================

# Using SparseArray.sparse matrix

sA = sparse([MI(1) MI(2); MI(3) MI(4)])
@test sA isa SparseMatrixCSC{MI, Int64}

sB = MI(sparse([1 2; 3 4]))
@test sB isa SparseMatrixCSC{MI, Int64}

sC = sparse([MI(1) 2; 3 4])
@test sC isa SparseMatrixCSC{MI, Int64}

sD = sparse([MI(1.0) 2; 3 4])
@test sD isa SparseMatrixCSC{MI, Int64}

sE = sparse([MI(1) 2.0; 3 4])
@test sE isa SparseMatrixCSC{MI, Int64}

sF = sparse([1 MI(2.0); 3 4])
@test sF isa SparseMatrixCSC{MI, Int64}

@test sA == sB == sC == sD == sE == sF

# Using SparseArray.sparse (I, J, D) vectors

spA = MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]))
@test size(spA.nzval,1) == 3
@test spA.nzval == [mp0; 2.0; 0.0]

spB = sparse(MP([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]))
@test size(spB.nzval,1) == 3
@test spB.nzval == [mp0; 2.0; 0.0]
@test spB == spA

spC = MP(spA)
@test size(spC.nzval,1) == 3
@test spC.nzval == [mp0; 2.0; 0.0]
@test spC == spA

# Using Max-Plus

spA = sparse(MP([-Inf 0; 0 -Inf]))
@test findnz(spA) == ([2, 1], [1, 2], MP[0.0, 0.0])
spB = MP(sparse([-Inf 0; 0 -Inf]))
@test findnz(spB) == ([1, 2], [1, 2], MP[mp0, mp0])
@test spA != spB
spC = sparse(MP([4 0; 7 -Inf]))
@test findnz(spC) == ([1, 2, 1], [1, 1, 2], MP[4, 7, 0])

# ==============================================================================
# Matrix ones, eye, zeros constructions
# ==============================================================================

# Max-Plus Matrix of ones

@test ones(MP, 2) isa Vector{MP}
@test ones(MP, 2) == [mp1; mp1]
@test ones(MP, 2,5) isa Matrix{MP}
@test ones(MP, 2,5) == [mp1 mp1 mp1 mp1 mp1; mp1 mp1 mp1 mp1 mp1]

# Max-Plus Identity matrix

@test eye(MP, 2) isa Matrix{MP}
@test eye(MP, 2) == [mp1 mp0; mp0 mp1]
@test eye(MP, 2,5) isa Matrix{MP}
@test eye(MP, 2,5) == [mp1 mp0 mp0 mp0 mp0; mp0 mp1 mp0 mp0 mp0]

# Max-Plus Matrix of zeros

@test spzeros(MP, 2) isa SparseVector{MP, Int64}
@test spzeros(MP, 2).nzval == MP([])
@test spzeros(MP, 2,3) isa SparseMatrixCSC{MP, Int64}
@test spzeros(MP, 2,3).nzval == MP([])

# Min-Plus Matrix of ones

@test ones(MI, 2) isa Vector{MI}
@test ones(MI, 2) == [mi1; mi1]
@test ones(MI, 2,5) isa Matrix{MI}
@test ones(MI, 2,5) == [mi1 mi1 mi1 mi1 mi1; mi1 mi1 mi1 mi1 mi1]

# Min-Plus Identity matrix

@test eye(MI, 2) isa Matrix{MI}
@test eye(MI, 2) == [mi1 mi0; mi0 mi1]
@test eye(MI, 2,5) isa Matrix{MI}
@test eye(MI, 2,5) == [mi1 mi0 mi0 mi0 mi0; mi0 mi1 mi0 mi0 mi0]

# Min-Plus Matrix of zeros

@test spzeros(MI, 2) isa SparseVector{MI, Int64}
@test spzeros(MI, 2).nzval == MI([])
@test spzeros(MI, 2,3) isa SparseMatrixCSC{MI, Int64}
@test spzeros(MI, 2,3).nzval == MI([])

# ==============================================================================
# Max-Plus matrix dense/sparse multiplication
# ==============================================================================

# Min-Plus TODO

# Max-Plus

A = [mp0 mp0 3; -2 mp0 mp0; mp0 mp0 0]
B = [mp0 2 mp0; mp0 mp0 0; -3 mp0 mp0]
D = [3 mp0 mp0; mp0 -2 mp0; mp0 mp0 0]
P = [mp0 mp0 0; 0 mp0 mp0; mp0 mp0 0]

@test A == (D * P)
@test A == (sparse(D) * P)
@test A == (D * sparse(P))
@test A == (sparse(D) * sparse(P))
@test sparse(A) == (sparse(D) * sparse(P))

# ==============================================================================
# Max-Plus matrix addition
# ==============================================================================

# Min-Plus TODO

# Max-Plus

A = MP([3.0 4; 5 6])
B = MP([1.0 2; 3 4])
C = MP([3.0 4; 5 6])
@test (A + B) == (B + A) == C

# ==============================================================================
# Matrix ones, eye, zeros operations
# ==============================================================================

# Max-Plus

I = eye(MP,4,4)
Z = spzeros(MP,4,4)
O = ones(MP,4,4)
A = MP(rand(4,4))

@test (A * I) == (I * A) == A
@test (A + Z) == (Z + A) == A
@test (A * Z) == (Z * A) == Z
@test (A + O) == (O + A) == A
@test (A * O) != (O * A) != A

# Min-Plus

I = eye(MI,4,4)
Z = spzeros(MI,4,4)
O = ones(MI,4,4)
A = MI(rand(4,4))

@test (A * I) == (I * A) == A
@test (A + Z) == (Z + A) == A
@test (A * Z) == (Z * A) == Z
@test (A + O) == (O + A) == O
@test (A * O) != (O * A) != A

# ==============================================================================
# Max-Plus element by element operations on matrices
# ==============================================================================

# Min-Plus TODO

# Max-Plus Dense matrix

A = MP([1.0 2.0; 3.0 4.0])
@test 2.0 .+ A == A .+ 2.0 == [MP(2.0) + MP(1.0) MP(2.0) + MP(2.0); MP(2.0) + MP(3.0) MP(2.0) + MP(4.0)] == MP([2.0 2.0; 3.0 4.0])
@test 2.0 .* A == A .* 2.0 == [MP(2.0) * MP(1.0) MP(2.0) * MP(2.0); MP(2.0) * MP(3.0) MP(2.0) * MP(4.0)] == MP([3.0 4.0; 5.0 6.0])

# Max-Plus Sparse matrix

sA = sparse(A)
@test 2.0 .+ sA == sA .+ 2.0 == [MP(2.0) + MP(1.0) MP(2.0) + MP(2.0); MP(2.0) + MP(3.0) MP(2.0) + MP(4.0)] == MP([2.0 2.0; 3.0 4.0])
@test 2.0 .* sA == sA .* 2.0 == [MP(2.0) * MP(1.0) MP(2.0) * MP(2.0); MP(2.0) * MP(3.0) MP(2.0) * MP(4.0)] == MP([3.0 4.0; 5.0 6.0])

# ==============================================================================
# Max-Plus type promotion/contamination on dense/sparse matrix
# ==============================================================================

# Dense matrix

A = [1 2; 3 4]
@test A isa Matrix{Int64}
@test MP(A) isa Matrix{MP}

B = [MP(1) MP(2); MP(3) MP(4)]
@test B isa Matrix{MP}

C = [MP(1) 2; 3 4]
@test C isa Matrix{MP}

D = [MP(1.0) 2; 3 4]
@test D isa Matrix{MP}

E = [MP(1) 2.0; 3 4]
@test E isa Matrix{MP}

F = [1 MP(2.0); 3 4]
@test F isa Matrix{MP}

# Sparse matrix

sA = sparse([1 2; 3 4])
@test sA isa SparseMatrixCSC{Int64, Int64}
@test MP(sA) isa SparseMatrixCSC{MP, Int64}

sB = sparse([MP(1) MP(2); MP(3) MP(4)])
@test sB isa SparseMatrixCSC{MP, Int64}

sC = sparse([MP(1) 2; 3 4])
@test sC isa SparseMatrixCSC{MP, Int64}

sD = sparse([MP(1.0) 2; 3 4])
@test sD isa SparseMatrixCSC{MP, Int64}

sE = sparse([MP(1) 2.0; 3 4])
@test sE isa SparseMatrixCSC{MP, Int64}

sF = sparse([1 MP(2.0); 3 4])
@test sF isa SparseMatrixCSC{MP, Int64}

# Dense/Sparse matrix comparaison

@test (A == sA) == (sA == A)
@test (B == sB) == (sB == B)
@test (C == sC) == (sC == C)
@test (D == sD) == (sD == D)
@test (E == sE) == (sE == E)
@test (F == sF) == (sF == F)

# ==============================================================================
# Max-Plus scalar to class algebra scalar conversion
# ==============================================================================

# Dense Matrix

@test plustimes(C) isa Matrix{Float64}
@test plustimes(D) isa Matrix{Float64}
@test plustimes(E) isa Matrix{Float64}
@test plustimes(F) isa Matrix{Float64}

# Sparse Matrix

@test plustimes(sC) isa SparseMatrixCSC{Float64, Int64}
@test plustimes(sD) isa SparseMatrixCSC{Float64, Int64}
@test plustimes(sE) isa SparseMatrixCSC{Float64, Int64}
@test plustimes(sF) isa SparseMatrixCSC{Float64, Int64}

# Max-Plus Sparse Matrix to Max-Plus Dense Matrix

@test full(sC) isa Matrix{MP}
@test full(sD) isa Matrix{MP}
@test full(sE) isa Matrix{MP}
@test full(sF) isa Matrix{MP}

@test full(sC) == C
@test full(sD) == D
@test full(sE) == E
@test full(sF) == F

# dense() is a anlias for full()

@test dense(sC) isa Matrix{MP}
@test dense(sD) isa Matrix{MP}
@test dense(sE) isa Matrix{MP}
@test dense(sF) isa Matrix{MP}

@test dense(sC) == C
@test dense(sD) == D
@test dense(sE) == E
@test dense(sF) == F

# Max-Plus Sparse Matrix to classic algebra Dense Matrix

@test array(sC) isa Matrix{Float64}
@test array(sD) isa Matrix{Float64}
@test array(sE) isa Matrix{Float64}
@test array(sF) isa Matrix{Float64}

# Max-Plus sparse array to Max-Plus dense array

Z = dense(spzeros(MP,2,2))
@test typeof(Z) == Matrix{MP}
@test Z == [mp0 mp0; mp0 mp0]

# Max-Plus sparse array to Max-Plus dense array

Z = full(spzeros(MP,2,2))
@test typeof(Z) == Matrix{MP}
@test Z == [mp0 mp0; mp0 mp0]

# Max-Plus sparse array to non Max-Plus dense array

Z = array(spzeros(MP,2,2))
@test typeof(Z) == Matrix{Float64}
@test Z == [-Inf -Inf; -Inf -Inf]

# ==============================================================================
# Bug with Julia with SparseMatrixCSC and operator== which confused zero() and 0.0
# ==============================================================================

# Min-Plus

A = MI(sparse([1, 2], [1, 2], [0.0, 0.0]))
B = spzeros(MI,2,2)
@test A.nzval == MI([0.0, 0.0])
@test (A == B) == false

AA = sparse([1, 2], [1, 2], [mi0, mi0])
BB = spzeros(MI,2,2)
@test AA.nzval == MI([Inf, Inf])
@test (AA == BB) == true

# Max-Plus

A = MP(sparse([1, 2], [1, 2], [0.0, 0.0]))
B = spzeros(MP,2,2)
@test A.nzval == MP([0.0, 0.0])
@test (A == B) == false

AA = sparse([1, 2], [1, 2], [mp0, mp0])
BB = spzeros(MP,2,2)
@test AA.nzval == MP([-Inf, -Inf])
@test (AA == BB) == true

# ==============================================================================
# Sparse Max-Plus element insertion
# ==============================================================================

# Note: ScicosLab will return isempty == false and length = 2
# Note: NSP, like Julia will return isempty == false and length = 4

O = spzeros(MP,2,2)
@test length(O) == 4
@test isempty(O) == false

# Zero elements are not inserted
O[1,1] = mp0
@test length(O) == 4
@test isempty(O) == false
@test O.nzval == []
@test isempty(nonzeros(O))

# Insert fake zero element
O[1,1] = MP(0.0)
@test O[1,1] == MP(0.0)
@test O.nzval == MP([0.0])
@test !isempty(nonzeros(O))

# Allow replacing a non-zero() element by a zero() element
O[1,1] = mp0
@test O[1,1] == MP(-Inf)
@test O.nzval == MP([-Inf])
@test !isempty(nonzeros(O))

# ==============================================================================
# Max-Plus scalar power operator
# ==============================================================================

A = MP([4 3; 7 -Inf])
@test A^0 == eye(MP,2,2)
@test A * A == A^2 == MP([10.0 7.0; 11.0 10.0])
@test A * A * A == A^3 == MP([14.0 13.0; 17.0 14.0])
@test_broken A^-1

# ==============================================================================
# Max-Plus trace
# ==============================================================================

A = [5 mp0 5; mp0 6 3; 11 12 11]

@test tr(A) == MP(11.0)
@test tr(MP([])) == mp0
@test tr(eye(MP,2,2)) == tr(eye(MP,2,2)) == 0.0
@test tr(spzeros(MP,2,2)) == tr(full(spzeros(MP,2,2))) == tr(spzeros(MP,2,2)) == mp0
@test tr(MP([1.0 2.0; 3.0 4.0])) == MP(1.0) + MP(4.0) == MP(4.0)
@test tr(sparse(MP([1.0 2.0; 3.0 4.0]))) == MP(4.0)

# ==============================================================================
# Max-Plus norm
# ==============================================================================

@test norm(MP([1.0 20 2; 30 400 4; 4 50 10])) == MP(400 - 1) == MP(399.0)
@test norm(MP([1 20 2; 30 400 4; 4 50 10])) == MP(400 - 1) == MP(399)
@test norm(sparse(MP([1 20 2; 30 400 4; 4 50 10]))) == MP(400 - 1) == MP(399)
@test norm([mp0 1; 10 mp0]) == MP(10.0 - -Inf) == mptop

# ==============================================================================
# Max-Plus inverse
# ==============================================================================

A = [mp0 1 mp0; 2 mp0 mp0; mp0 mp0 3]
@test inv(A) == A^-1 == [mp0 -2 mp0; -1 mp0 mp0; mp0 mp0 -3]
@test A * inv(A) == inv(A) * A == eye(MP,3,3)
@test typeof(inv(A)) == typeof(A^-1) == Matrix{MP}

A = [mp0 1 mp0; 2 mp0 mp0]
@test inv(A) == A^-1 == [mp0 -2; -1 mp0; mp0 mp0]
@test A * inv(A) == eye(MP,2,2)
@test inv(A) * A == [0 mp0 mp0; mp0 0 mp0; mp0 mp0 mp0]
@test A * inv(A) != inv(A) * A

#FIXME @test_throws ErrorException("The matrix cannot be inversed") inv(MP([1 2; 3 4]))

# ==============================================================================
# Matrix residuation
# ==============================================================================

A = [mp0 1 mp0; 2 mp0 mp0; mp0 mp0 3]
B = [3 mp0 mp0; mp0 mp0 4; mp0 5 mp0]
x = A \ B
@test x == MP([mp0 mp0 2; 2 mp0 mp0; mp0 2 mp0])
@test A * x == B
@test (A \ A) == (A / A) == eye(MP,3,3)
@test (B \ B) == (B / B) == eye(MP,3,3)

###

A = MP([3.0 4; 5 6])
B = MP([0.0 -2; 2 0])
x = A \ B
@test x == MP([-3 -5; -4 -6])
@test A * x == B

@test (A \ A) == MP([0 1; -1 0])
@test (A / A) == (B \ B) == (B / B) == MP([0 -2; 2 0])

# ==============================================================================
# Max-Plus star
# ==============================================================================

# Scalars

@test star(MP(2)) == mptop
@test star(MP(1.0)) == mptop
@test star(MP(-1.0)) == mp1
@test star(mp0) == mp1
@test star(mp1) == mp1
@test star(mptop) == mptop

# Matrices

@test_throws ErrorException("Matrix shall be squared") star(MP([]))
@test_throws ErrorException("Matrix shall be squared") star(eye(MP,3,2))
@test star(eye(MP,2,2)) == eye(MP,2,2)
@test star(full(spzeros(MP,2,2))) == eye(MP,2,2)

#

@test star(MP([1 2; 3 4])) == [mptop mptop; mptop mptop]
A = MP([-3 -2; -1 0]); B = star(A)
@test B == eye(MP,2,2) + A
@test B == B + A^2
@test B == B + A^3

#

@test star(MP([-1 2; mp0 -3])) == MP([0 2; -Inf 0])
A = [mp0 2 3; -2 -10 -1; -5 -2 mp1]
@test star(A) == MP([0 2 3; -2 0 1; -4 -2 0])
@test star(star(A)) == star(A)
@test star(A) == (A^0 + A)^2
@test (A^0 + A)^2 == (A^0 + A)^3

# FIXME KO: Mixing +inf and -inf

# @test mpstar(MP([2 3; mp0 -1])) == [mi0 mi0; mp0 mp1]

# Random large matrix

A = MP(rand(64,64))
@test star(A) == fill(mptop, 64,64)

# FIXME KO
B = (((ones(1, size(A,1)) * A * ones(size(A,2), 1))[1,1])^-1) * A
@test_broken maximum(plustimes(B)) == 0.0

# ==============================================================================
# Max-Plus plus
# ==============================================================================

# Scalars

@test plus(MP(2)) == mptop
@test plus(MP(1.0)) == mptop
@test_broken plus(MP(-1.0)) == MP(-1)
@test_broken plus(mp0) == mp0
@test plus(mp1) == mp1
@test plus(mptop) == mptop

# Matrices

@test_throws ErrorException("Matrix shall be squared") plus(MP([]))
@test_throws ErrorException("Matrix shall be squared") plus(eye(MP,3,2))
A = [mp0 2 3; -2 -10 -1; -5 -2 mp1]
B = plus(A)
@test B == MP([0 2 3; -2 0 1; -4 -2 0])
@test B == A * star(A)

# FIXME donne le bon resultat dans REPL
A[2,1] = MP(-10)
@test_broken plus(A) == MP([-2 2 3; -6 -3 -1; -5 -2 0])

# What happens if a circuit has strictly positive weight ?
A[3,1] = 6
@test_broken plus(A) == fill(mi0, 2,3) # FIXME

# ==============================================================================
# TODO Max-Plus a star b
# ==============================================================================

# ==============================================================================
# Max-Plus Howard
# ==============================================================================

S = sparse(MP([1 2; 3 4]))
λ,v = howard(S)
@test (λ,v) == (MP[4, 4], MP[2, 4])
@test (S * v) == (λ[1] * v)

S = sparse([mp0 2 mp0; mp1 mp0 mp0; mp0 mp0 2])
λ,v = howard(S)
@test (λ,v) == (MP[1, 1, 2], MP[1, 0, 2])
B = S / Matrix(Diagonal(λ))
@test B == [mp0 1 mp0; -1 mp0 mp0; mp0 mp0 0]
@test B * v == v

S = sparse([1 1; mp0 2])
λ,v = howard(S)
@test (λ,v) == (MP[2, 2], MP[1, 2])
@test (S * v) == (λ[1] * v)
@test (S * [0; mp0]) == (MP(1) * [0; mp0])

S = sparse([2 1; mp0 mp1])
λ,v = howard(S)
@test (λ,v) == (MP[2, 0], MP[2, 0])
B = S / Matrix(Diagonal(λ))
@test B == [0 1; mp0 0]
@test B * v == v

# ==============================================================================
# Display
# ==============================================================================

using Suppressor

# Min-Plus TODO

# Max-Plus

mp_change_display(0)
result = @capture_out show(stdout, mp0)
@test result == "-Inf"
result = @capture_out show(stdout, mp1)
@test result == "0.0"
result = @capture_out show(stdout, MP(4.0))
@test result == "4.0"
result = @capture_out show(stdout, zero(MP))
@test result == "-Inf"
result = @capture_out show(stdout, one(MP))
@test result == "0.0"

mp_change_display(1)
result = @capture_out show(stdout, mp0)
@test result == "."
result = @capture_out show(stdout, mp1)
@test result == "0"
result = @capture_out show(stdout, MP(4.0))
@test result == "4"
result = @capture_out show(stdout, zero(MP))
@test result == "."
result = @capture_out show(stdout, one(MP))
@test result == "0"

mp_change_display(2)
result = @capture_out show(stdout, mp0)
@test result == "."
result = @capture_out show(stdout, mp1)
@test result == "e"
result = @capture_out show(stdout, MP(4.0))
@test result == "4"
result = @capture_out show(stdout, zero(MP))
@test result == "."
result = @capture_out show(stdout, one(MP))
@test result == "e"

mp_change_display(3)
result = @capture_out show(stdout, mp0)
@test result == "ε"
result = @capture_out show(stdout, mp1)
@test result == "0"
result = @capture_out show(stdout, MP(4))
@test result == "4"
result = @capture_out show(stdout, zero(MP))
@test result == "ε"
result = @capture_out show(stdout, one(MP))
@test result == "0"

mp_change_display(4)
result = @capture_out show(stdout, mp0)
@test result == "ε"
result = @capture_out show(stdout, mp1)
@test result == "e"
result = @capture_out show(stdout, MP(4.0))
@test result == "4"
result = @capture_out show(stdout, zero(MP))
@test result == "ε"
result = @capture_out show(stdout, one(MP))
@test result == "e"

A = MP([4.5 0.0; 7.0 -Inf])
result = @capture_out show(stdout, A)
@test result == "MP[4.5 e; 7 ε]"
result = @capture_out LaTeX(stdout, A)
@test result == "\\left[\n\\begin{array}{*{20}c}\n4.5 & e \\\\\n7 & \\varepsilon \\\\\n\\end{array}\n\\right]\n"
mp_change_display(0)
result = @capture_out LaTeX(stdout, A)
@test result == "\\left[\n\\begin{array}{*{20}c}\n4.5 & 0 \\\\\n7 & -\\infty \\\\\n\\end{array}\n\\right]\n"


text:
# Ideas from MathML.jl

"""
Map XML tag names to parser functions.

$(TYPEDEF)
"""
const tagmap = Dict{String,Function}(
    "and" => parse_and,
    "arbitraryoperator" => parse_arbitraryoperator,
    "arbitrarysort" => parse_arbitrarysort,
    "arc"  => parse_arc,
    "bool" => parse_bool,
    "booleanconstant" => parse_booleanconstant,
    "condition" => parse_condition,
    "declaration" => parse_declaration,
    "declarations"  => parse_declarations,
    "equality" => parse_equality,
    "graphics" => parse_graphics,
    "hlinitialMarking" => parse_hlinitialMarking,
    "hlinscription" => parse_hlinscription,
    "imply" => parse_imply,
    "inequality" => parse_inequality,
    "initialMarking" => parse_initialMarking,
    "inscription" => parse_inscription,
    "label" => parse_label,
    "mulitsetsort" => parse_mulitsetsort,
    "name" => parse_name,
    "namedoperator" => parse_namedoperator,
    #"net" => parse_net, #! Method should not be forun by parse_node.
    "not" => parse_not,
    "or" => parse_or,
    "page"  => parse_page,
    "place"  => parse_place,
    #"pnml" => parse_pnml, #! Method should not be forun by parse_node.
    "productsort" => parse_productsort,
    "referencePlace"  => parse_refPlace,
    "referenceTransition"  => parse_refTransition,
    "sort" => parse_sort,
    "structure" => parse_structure,
    "term" => parse_term, #TODO is this valid?
    "text" => parse_text,
    "tokengraphics"  => parse_tokengraphics,
    "tokenposition" => parse_tokenposition,
    "toolspecific"  => parse_toolspecific,
    "transition"  => parse_transition,
    "tuple" => parse_tuple,
    "type" => parse_type,
    "unparsed" => parse_unparsed,
    "useroperator" => parse_useroperator,
    "usersort" => parse_usersort,
    "variable" => parse_variable,
    "variabledecl" => parse_variabledecl,

    # "add" => unclaimed_label,
    # "addition" => unclaimed_label,
    # "all" => unclaimed_label,
    # "anyName" => unclaimed_label,
    # "arctype" => unclaimed_label,
    # "attribute" => unclaimed_label,
    # "bool" => unclaimed_label,
    # "cardinality" => unclaimed_label,
    # "cardinalityof" => unclaimed_label,
    # "choice" => unclaimed_label,
    # "cyclicenumeration" => unclaimed_label,
    # "data" => unclaimed_label,
    # "def" => unclaimed_label,
    # "define" => unclaimed_label,
    # "dimension" => unclaimed_label,
    # "div" => unclaimed_label,
    # "dot" => unclaimed_label,
    # "dotconstant" => unclaimed_label,
    # "element" => unclaimed_label,
    # "empty" => unclaimed_label,
    # "emptylist" => unclaimed_label,
    # "except" => unclaimed_label,
    # "externalRef" => unclaimed_label,
    # "feconstant" => unclaimed_label,
    # "fill" => unclaimed_label,
    # "finiteenumeration" => unclaimed_label,
    # "finiteintrange" => unclaimed_label,
    # "finiteintrangeconstant" => unclaimed_label,
    # "font" => unclaimed_label,
    # "geq" => unclaimed_label,
    # "geqs" => unclaimed_label,
    # "grammar" => unclaimed_label,
    # "greaterthan" => unclaimed_label,
    # "greaterthanorequal" => unclaimed_label,
    # "group" => unclaimed_label,
    # "gt" => unclaimed_label,
    # "gtp" => unclaimed_label,
    # "gts" => unclaimed_label,
    # "include" => unclaimed_label,
    # "input" => unclaimed_label,
    # "integer" => unclaimed_label,
    # "interleave" => unclaimed_label,
    # "leq" => unclaimed_label,
    # "leqs" => unclaimed_label,
    # "lessthan" => unclaimed_label,
    # "lessthanorequal" => unclaimed_label,
    # "line" => unclaimed_label,
    # "list" => unclaimed_label,
    # "listappend" => unclaimed_label,
    # "listconcatenation" => unclaimed_label,
    # "listlength" => unclaimed_label,
    # "lt" => unclaimed_label,
    # "ltp" => unclaimed_label,
    # "lts" => unclaimed_label,
    # "makelist" => unclaimed_label,
    # "memberatindex" => unclaimed_label,
    # "mixed" => unclaimed_label,
    # "mod" => unclaimed_label,
    # "mult" => unclaimed_label,
    # "multisetsort" => unclaimed_label,
    # "namedsort" => unclaimed_label,
    # "natural" => unclaimed_label,
    # "notAllowed" => unclaimed_label,
    # "nsName" => unclaimed_label,
    # "numberconstant" => unclaimed_label,
    # "numberof" => unclaimed_label,
    # "offset" => unclaimed_label,
    # "oneOrMore" => unclaimed_label,
    # "optional" => unclaimed_label,
    # "output" => unclaimed_label,
    # "param" => unclaimed_label,
    # "parameter" => unclaimed_label,
    # "parentRef" => unclaimed_label,
    # "partition" => unclaimed_label,
    # "partitionelement" => unclaimed_label,
    # "partitionelementof" => unclaimed_label,
    # "position" => unclaimed_label,
    # "positive" => unclaimed_label,
    # "predecessor" => unclaimed_label,
    # "ref" => unclaimed_label,
    # "scalarproduct" => unclaimed_label,
    # "start" => unclaimed_label,
    # "string" => unclaimed_label,
    # "stringappend" => unclaimed_label,
    # "stringconcatenation" => unclaimed_label,
    # "stringconstant" => unclaimed_label,
    # "stringlength" => unclaimed_label,
    # "sublist"  => unclaimed_label,
    # "substring" => unclaimed_label,
    # "subterm" => unclaimed_label,
    # "subtract" => unclaimed_label,
    # "subtraction" => unclaimed_label,
    # "successor" => unclaimed_label,
    # "value" => unclaimed_label,
    # "zeroOrMore" => unclaimed_label,
)

ext: jl


text:
# Define new Python classes from Julia.
#############################################################################

import MacroTools: @capture, splitdef, combinedef

######################################################################
# def_py_class definition: this is the core non-macro interface for creating
# a Python class from a Julia type.

"""
    def_py_class(type_name::AbstractString, methods::Vector;
                 base_classes=[], getsets::Vector=[])

`def_py_class` creates a Python class whose methods are implemented in Julia.
`@pydef` macros expand into a call to this function.

Arguments
---------
- `methods`: a vector of tuples `(py_name::String, jl_fun::Function)`
   `py_name` will be a method of the Python class, which will call `jl_fun`
- `base_classes`: the Python base classes to inherit from.

Return value: the created class (`::PyTypeObject`)
"""
function def_py_class(type_name::AbstractString, methods::Vector;
                      base_classes=[], getsets::Vector=[])
    # Only create new-style classes
    base_classes = union(base_classes, [pybuiltin("object")])
    methods = Dict(py_name => jlfun2pyfun(jl_fun::Function)
                   for (py_name::Symbol, jl_fun) in methods)
    getter_setters = Dict(py_name => pyproperty(jlfun2pyfun(getter),
                                                jlfun2pyfun(setter))
                          for (py_name::Symbol, getter, setter) in getsets)
    return pybuiltin("type")(type_name, tuple(base_classes...),
                             merge(methods, getter_setters))
end

######################################################################
# @pydef macro

# Helper for `parse_pydef`
# Returns (class_name::Symbol, base_classes, lines)
# where there's one `line` per method definition
function parse_pydef_toplevel(expr)
    if @capture(expr, begin mutable struct class_name_ <: base_classes_expr_
                    lines__
                end end)
        if !@capture(base_classes_expr, (base_classes__,))
            base_classes = (base_classes_expr,)
        end
    else
        @assert(@capture(expr, mutable struct class_name_
                    lines__
            end), "Malformed @pydef expression")

        base_classes = []
    end
    return class_name::Symbol, base_classes, lines
end

function parse_pydef(expr)
    class_name, base_classes, lines = parse_pydef_toplevel(expr)
    # Now we parse every method definition / getter / setter
    function_defs = Expr[] # vector of :(function ...) expressions
    methods = Tuple{Any, Symbol}[] # (py_name, jl_method_name)
    getter_dict = Dict{Any, Symbol}() # python_var => jl_getter_name
    setter_dict = Dict{Any, Symbol}()
    method_syms = Dict{Any, Symbol}() # see below
    for line in lines
        if !isa(line, LineNumberNode) && line.head != :line # need to skip those
            def_dict = splitdef(line)
            py_f = def_dict[:name]
            # The dictionary of the new Julia-side definition.
            jl_def_dict = copy(def_dict)
            if isa(py_f, Symbol)
                # Method definition
                # We save the gensym to support multiple dispatch
                #    readlines(io) = ...
                #    readlines(io, nlines) = ...
                # otherwise the first and second `readlines` get different
                # gensyms, and one of the two gets shadowed by the other.
                jl_def_dict[:name] = get!(method_syms, py_f, gensym(py_f))
                if py_f == :__init__
                    # __init__ must return `nothing` in Python. This is
                    # achieved by default in Python, but not so in Julia, so we
                    # special-case it for convenience.
                    jl_def_dict[:body] = :(begin $(def_dict[:body]); nothing end)
                end
                push!(methods, (py_f, jl_def_dict[:name]))
            elseif @capture(py_f, attribute_.access_)
                # Accessor (.get/.set) definition
                if access == :get
                    dict = getter_dict
                elseif access == :set!
                    dict = setter_dict
                else
                    error("$access is not a valid accessor; must be either get or set!")
                end
                dict[attribute] = jl_def_dict[:name] = gensym(Symbol(attribute,:_,access))
            else
                error("Malformed line: $line")
            end
            push!(function_defs, combinedef(jl_def_dict))
        end
    end
    @assert(isempty(setdiff(keys(setter_dict), keys(getter_dict))),
            "All .set attributes must have a .get")
    class_name, base_classes, methods, getter_dict, setter_dict, function_defs
end

"""
`@pydef` creates a Python class whose methods are implemented in Julia.
For instance,

    @pyimport numpy.polynomial as P
    @pydef type Doubler <: P.Polynomial
        __init__(self, x=10) = (self[:x] = x)
        my_method(self, arg1::Number) = arg1 + 20
        x2.get(self) = self[:x] * 2
        x2.set!(self, new_val) = (self[:x] = new_val / 2)
    end
    Doubler()[:x2]

is essentially equivalent to the following Python code:

    class Doubler(numpy.polynomial.Polynomial):
        def __init__(self, x=10):
            self.x = x
        def my_method(self, arg1): return arg1 + 20
        @property
        def x2(self): return self.x * 2
        @x2.setter
        def x2(self, new_val):
            self.x = new_val / 2
    Doubler().x2

The method arguments and return values are automatically converted between Julia and Python. All Python
special methods are supported (`__len__`, `__add__`, etc.)

`@pydef` allows for multiple inheritance of Python types:

    @pydef type SomeType <: (BaseClass1, BaseClass2)
        ...
    end

Multiple dispatch works, too:

    x2.set!(self, new_x::Int) = ...
    x2.set!(self, new_x::Float64) = ...
"""
macro pydef(class_expr)
    class_name, _, _ = parse_pydef_toplevel(class_expr)
    esc(:(const $class_name = $PyCall.@pydef_object($class_expr)))
end

"""
`@pydef_object` is like `@pydef`, but it returns the
metaclass as a `PyObject` instead of binding it to the class name.
It's side-effect-free, except for the definition of the class methods.
"""
macro pydef_object(class_expr)
    class_name, base_classes, methods_, getter_dict, setter_dict, function_defs=
        parse_pydef(class_expr)
    methods = [:($(Expr(:quote, py_name::Symbol)), $(esc(jl_fun::Symbol)))
               for (py_name, jl_fun) in methods_]
    getsets = [:($(Expr(:quote, attribute)),
                 $(esc(getter)),
                 $(esc(get(setter_dict, attribute, nothing))))
               for (attribute, getter) in getter_dict]
    :(begin
        $(map(esc, function_defs)...)
        # This line doesn't have any side-effect, it just returns the Python
        # (meta-)class, as a PyObject
        def_py_class($(string(class_name)), [$(methods...)];
                     base_classes=[$(map(esc, base_classes)...)],
                     getsets=[$(getsets...)])
    end)
end


text:
####
# FB Oracle tests
####

facts("suppFcnTest FB") do
	mfs     = [.1, .2]
	mbs     = [-.05, -.25]
	sigfs   = [1., 2.]
	sigbs   = [2., .5]
	w       = FBOracle(mfs, mbs, sigfs, sigbs, .1)

	zstar, ustar = suppFcn([1, 1], w, :Min)
	@fact zstar -->  roughly(-4.724022297688992)
	@fact ustar[1] --> roughly(-4.213785691942581)
	@fact ustar[2] --> roughly(-0.5102366057464114)

	zstar, ustar = suppFcn([1, 1], w, :Max)
	@fact zstar --> roughly(5.098525912188081)
	@fact ustar[1] --> roughly(1.0597051824376162)
	@fact ustar[2] --> roughly(4.038820729750465)

	zstar, ustar = suppFcn([-1,  1], w, :Min)
	@fact zstar --> roughly(-2.7492629560940403)
	@fact ustar[1] --> roughly(2.0194103648752324)
	@fact ustar[2] --> roughly(-0.7298525912188081)

	zstar, ustar = suppFcn([-1,  1], w, :Max)
	@fact zstar --> roughly(6.319708517540586)
	@fact ustar[1] --> roughly(-3.0848542587702927)
	@fact ustar[2] --> roughly(3.234854258770293)
end

facts("portTest for FB") do
	srand(8675309); data = randn(500, 2)
	w = FBOracle(data, .1, .2)
	portTest(w, -1.791121672173153, [0.5840372921578592, 0.41596270784214084],TOL=1e-2)
end


text:
# -----------------------------------------------------------------------------
##### CUBOID #####-------------------------------------------------------------
function subregioncuboid(dataobject::ClumpDataType;
        xrange::Array{<:Any,1}=[missing, missing],
        yrange::Array{<:Any,1}=[missing, missing],
        zrange::Array{<:Any,1}=[missing, missing],
        center::Array{<:Any,1}=[0., 0., 0.],
        range_unit::Symbol=:standard,
        inverse::Bool=false,
        verbose::Bool=verbose_mode)

    printtime("", verbose)

       boxlen = dataobject.boxlen

       # convert given ranges and print overview on screen
       ranges = prepranges(dataobject.info,range_unit, verbose, xrange, yrange, zrange, center)

       xmin, xmax, ymin, ymax, zmin, zmax = ranges

       #if !(xrange == [dataobject.ranges[1], dataobject.ranges[2]] &&
       #  yrange == [dataobject.ranges[3], dataobject.ranges[4]] &&
       #  zrange == [dataobject.ranges[5], dataobject.ranges[6]])
       if !(xrange[1] === missing && xrange[2] === missing &&
            yrange[1] === missing && yrange[2] === missing &&
            zrange[1] === missing && zrange[2] === missing)

          if inverse == false
              sub_data = filter(p->   p.peak_x >=  xmin * boxlen  &&
                                      p.peak_x <=  xmax * boxlen  &&
                                      p.peak_y >=  ymin * boxlen  &&
                                      p.peak_y <=  ymax * boxlen  &&
                                      p.peak_z >=  zmin * boxlen  &&
                                      p.peak_z <=  zmax * boxlen, dataobject.data)
          elseif inverse == true
              sub_data = filter(p->   (p.peak_x <  xmin * boxlen  ||
                                      p.peak_x >  xmax * boxlen)  ||
                                      (p.peak_y <  ymin * boxlen  ||
                                      p.peak_y >  ymax * boxlen)  ||
                                      (p.peak_z <  zmin * boxlen  ||
                                      p.peak_z >  zmax * boxlen), dataobject.data)
              ranges = dataobject.ranges
          end

          printtablememory(sub_data, verbose)

          clumpdata = ClumpDataType()
          clumpdata.data = sub_data
          clumpdata.info = dataobject.info
          clumpdata.boxlen = dataobject.boxlen
          clumpdata.ranges = ranges
          clumpdata.selected_clumpvars = dataobject.selected_clumpvars
          clumpdata.used_descriptors = dataobject.used_descriptors
          clumpdata.scale = dataobject.scale
          return clumpdata
      else
          return dataobject
      #    println("[Mera]: Nothing to do! Given ranges match data ranges!")
      #    println()
      end
end

# -----------------------------------------------------------------------------
##### CYLINDER #####-----------------------------------------------------------
function subregioncylinder(dataobject::ClumpDataType;
                            radius::Real=0.,
                            height::Real=0.,
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            direction::Symbol=:z,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    if radius == 0. || height == 0. || in(0., center)
        error("[Mera]: given radius, height or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_shift, height_shift = prepranges(dataobject.info, center, radius, height, range_unit, verbose)

    if inverse == false
        sub_data = filter(p-> sqrt( (p.peak_x -  cx_shift*boxlen)^2 +
                                    (p.peak_y -  cy_shift*boxlen )^2)
                                    <= ( radius_shift*boxlen )  &&
                            abs(p.peak_z - cz_shift*boxlen) <= ( height_shift*boxlen),
                                dataobject.data)
    elseif inverse == true
        sub_data = filter(p-> sqrt( (p.peak_x -  cx_shift*boxlen)^2 +
                                    (p.peak_y -  cy_shift*boxlen )^2)
                                    > ( radius_shift*boxlen )  ||
                            abs(p.peak_z - cz_shift*boxlen) > ( height_shift*boxlen),
                                dataobject.data)
        ranges = dataobject.ranges
    end

    printtablememory(sub_data, verbose)

    clumpdata = ClumpDataType()
    clumpdata.data = sub_data
    clumpdata.info = dataobject.info
    clumpdata.boxlen = dataobject.boxlen
    clumpdata.ranges = ranges
    clumpdata.selected_clumpvars = dataobject.selected_clumpvars
    clumpdata.used_descriptors = dataobject.used_descriptors
    clumpdata.scale = dataobject.scale
    return clumpdata
end

# -----------------------------------------------------------------------------
##### SPHERE #####-------------------------------------------------------------
function subregionsphere(dataobject::ClumpDataType;
                            radius::Real=0.,
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    if radius == 0. || in(0., center)
        error("[Mera]: given radius or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale

    # convert given ranges and print overview on screen
    height = 0.
    ranges, cx_shift, cy_shift, cz_shift, radius_shift = prepranges(dataobject.info, center, radius, height, range_unit, verbose)

    if inverse == false
        sub_data = filter(p-> sqrt( (p.peak_x -  cx_shift*boxlen)^2 +
                                    (p.peak_y -  cy_shift*boxlen )^2+
                                    (p.peak_z -  cz_shift*boxlen)^2 )
                                    <= ( radius_shift*boxlen ),
                                dataobject.data)
    elseif inverse == true
        sub_data = filter(p-> sqrt( (p.peak_x -  cx_shift*boxlen)^2 +
                                    (p.peak_y -  cy_shift*boxlen )^2+
                                    (p.peak_z -  cz_shift*boxlen)^2 )
                                    > ( radius_shift*boxlen ),
                                dataobject.data)
        ranges = dataobject.ranges
    end

    printtablememory(sub_data, verbose)

    clumpdata = ClumpDataType()
    clumpdata.data = sub_data
    clumpdata.info = dataobject.info
    clumpdata.boxlen = dataobject.boxlen
    clumpdata.ranges = ranges
    clumpdata.selected_clumpvars = dataobject.selected_clumpvars
    clumpdata.used_descriptors = dataobject.used_descriptors
    clumpdata.scale = dataobject.scale
    return clumpdata
end


text:
#########
# Trial #
#########

mutable struct Trial
    params::Parameters
    times::Vector{Float64}
    gctimes::Vector{Float64}
    memory::Vector{Int}
    allocs::Vector{Int}
end
StructTypes.StructType(::Type{Trial}) = StructTypes.Mutable()

Trial(params::Parameters) = Trial(params, Float64[], Float64[], Int[], Int[])
Trial() = Trial(Parameters())

function Base.:(==)(a::Trial, b::Trial)
    return a.params == b.params &&
           a.times == b.times &&
           a.gctimes == b.gctimes &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::Trial) = Trial(copy(t.params), copy(t.times), copy(t.gctimes), copy(t.memory), copy(t.allocs))

function Base.push!(t::Trial, time, gctime, memory, allocs)
    push!(t.times, time)
    push!(t.gctimes, gctime)
    push!(t.memory, memory)
    push!(t.allocs, allocs)
    return t
end

function Base.deleteat!(t::Trial, i)
    deleteat!(t.times, i)
    deleteat!(t.gctimes, i)
    deleteat!(t.memory, i)
    deleteat!(t.allocs, i)
    return t
end

Base.length(t::Trial) = length(t.times)
Base.getindex(t::Trial, i::Number) = push!(Trial(t.params), t.times[i], t.gctimes[i], t.memory[i], t.allocs[i])
Base.getindex(t::Trial, i) = Trial(t.params, t.times[i], t.gctimes[i], t.memory[i], t.allocs[i])
Base.lastindex(t::Trial) = length(t)

function Base.sort!(t::Trial)
    inds = sortperm(t.times)
    t.times = t.times[inds]
    t.gctimes = t.gctimes[inds]
    t.memory = t.memory[inds]
    t.allocs = t.allocs[inds]
    return t
end

Base.sort(t::Trial) = sort!(copy(t))

Base.time(t::Trial) = time(minimum(t))
gctime(t::Trial) = gctime(minimum(t))
memory(t::Trial) = memory(minimum(t))
allocs(t::Trial) = allocs(minimum(t))
params(t::Trial) = t.params

# returns the index of the first outlier in `values`, if any outliers are detected.
# `values` is assumed to be sorted from least to greatest, and assumed to be right-skewed.
function skewcutoff(values)
    current_values = copy(values)
    while mean(current_values) > median(current_values)
        deleteat!(current_values, length(current_values))
    end
    return length(current_values) + 1
end

skewcutoff(t::Trial) = skewcutoff(t.times)

function rmskew!(t::Trial)
    sort!(t)
    i = skewcutoff(t)
    i <= length(t) && deleteat!(t, i:length(t))
    return t
end

function rmskew(t::Trial)
    st = sort(t)
    return st[1:(skewcutoff(st) - 1)]
end

trim(t::Trial, percentage = 0.1) = t[1:max(1, floor(Int, length(t) - (length(t) * percentage)))]

#################
# TrialEstimate #
#################

mutable struct TrialEstimate
    params::Parameters
    time::Float64
    gctime::Float64
    memory::Float64
    allocs::Float64
end
StructTypes.StructType(::Type{TrialEstimate}) = StructTypes.Mutable()

function TrialEstimate(trial::Trial, t, gct, mem, alloc)
    return TrialEstimate(params(trial), t, gct, mem, alloc)
end

function Base.:(==)(a::TrialEstimate, b::TrialEstimate)
    return a.params == b.params &&
           a.time == b.time &&
           a.gctime == b.gctime &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::TrialEstimate) = TrialEstimate(copy(t.params), t.time, t.gctime, t.memory, t.allocs)

function Base.minimum(trial::Trial)
    i = argmin(trial.times)
    return TrialEstimate(trial, trial.times[i], trial.gctimes[i], trial.memory[i], trial.allocs[i])
end

function Base.maximum(trial::Trial)
    i = argmax(trial.times)
    return TrialEstimate(trial, trial.times[i], trial.gctimes[i], trial.memory[i], trial.allocs[i])
end

Statistics.median(trial::Trial) = TrialEstimate(trial, median(trial.times), median(trial.gctimes), median(trial.memory), median(trial.allocs))
Statistics.mean(trial::Trial) = TrialEstimate(trial, mean(trial.times), mean(trial.gctimes), mean(trial.memory), mean(trial.allocs))
Statistics.std(trial::Trial) = TrialEstimate(trial, std(trial.times), std(trial.gctimes), std(Float64.(trial.memory)), std(Float64.(trial.allocs)))

Base.isless(a::TrialEstimate, b::TrialEstimate) = isless(time(a), time(b))

Base.time(t::TrialEstimate) = t.time
gctime(t::TrialEstimate) = t.gctime
memory(t::TrialEstimate) = t.memory
allocs(t::TrialEstimate) = t.allocs
iallocs(t::TrialEstimate) = round(Int, t.allocs)
params(t::TrialEstimate) = t.params

##############
# TrialRatio #
##############

mutable struct TrialRatio
    params::Parameters
    time::Float64
    gctime::Float64
    memory::Float64
    allocs::Float64
end
StructTypes.StructType(::Type{TrialRatio}) = StructTypes.Mutable()

function Base.:(==)(a::TrialRatio, b::TrialRatio)
    return a.params == b.params &&
           a.time == b.time &&
           a.gctime == b.gctime &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::TrialRatio) = TrialRatio(copy(t.params), t.time, t.gctime, t.memory, t.allocs)

Base.time(t::TrialRatio) = t.time
gctime(t::TrialRatio) = t.gctime
memory(t::TrialRatio) = t.memory
allocs(t::TrialRatio) = t.allocs
params(t::TrialRatio) = t.params

function ratio(a::Real, b::Real)
    if a == b # so that ratio(0.0, 0.0) returns 1.0
        return one(Float64)
    end
    return Float64(a / b)
end

function ratio(a::TrialEstimate, b::TrialEstimate)
    ttol = max(params(a).time_tolerance, params(b).time_tolerance)
    mtol = max(params(a).memory_tolerance, params(b).memory_tolerance)
    p = Parameters(params(a); time_tolerance = ttol, memory_tolerance = mtol)
    return TrialRatio(p, ratio(time(a), time(b)), ratio(gctime(a), gctime(b)),
                      ratio(memory(a), memory(b)), ratio(allocs(a), allocs(b)))
end

gcratio(t::TrialEstimate) =  ratio(gctime(t), time(t))

##################
# TrialJudgement #
##################

struct TrialJudgement
    ratio::TrialRatio
    time::Symbol
    memory::Symbol
end
StructTypes.StructType(::Type{TrialJudgement}) = StructTypes.Struct()

function TrialJudgement(r::TrialRatio)
    ttol = params(r).time_tolerance
    mtol = params(r).memory_tolerance
    return TrialJudgement(r, judge(time(r), ttol), judge(memory(r), mtol))
end

function Base.:(==)(a::TrialJudgement, b::TrialJudgement)
    return a.ratio == b.ratio &&
           a.time == b.time &&
           a.memory == b.memory
end

Base.copy(t::TrialJudgement) = TrialJudgement(copy(t.params), t.time, t.memory)

Base.time(t::TrialJudgement) = t.time
memory(t::TrialJudgement) = t.memory
ratio(t::TrialJudgement) = t.ratio
params(t::TrialJudgement) = params(ratio(t))

judge(a::TrialEstimate, b::TrialEstimate; kwargs...) = judge(ratio(a, b); kwargs...)

function judge(r::TrialRatio; kwargs...)
    newr = copy(r)
    newr.params = Parameters(params(r); kwargs...)
    return TrialJudgement(newr)
end

function judge(ratio::Real, tolerance::Float64)
    if isnan(ratio) || (ratio - tolerance) > 1.0
        return :regression
    elseif (ratio + tolerance) < 1.0
        return :improvement
    else
        return :invariant
    end
end

isimprovement(f, t::TrialJudgement) = f(t) == :improvement
isimprovement(t::TrialJudgement) = isimprovement(time, t) || isimprovement(memory, t)

isregression(f, t::TrialJudgement) = f(t) == :regression
isregression(t::TrialJudgement) = isregression(time, t) || isregression(memory, t)

isinvariant(f, t::TrialJudgement) = f(t) == :invariant
isinvariant(t::TrialJudgement) = isinvariant(time, t) && isinvariant(memory, t)

const colormap = (
    regression = :red,
    improvement = :green,
    invariant = :normal,
)

printtimejudge(io, t::TrialJudgement) =
    printstyled(io, time(t); color=colormap[time(t)])
printmemoryjudge(io, t::TrialJudgement) =
    printstyled(io, memory(t); color=colormap[memory(t)])

###################
# Pretty Printing #
###################

prettypercent(p) = string(@sprintf("%.2f", p * 100), "%")

function prettydiff(p)
    diff = p - 1.0
    return string(diff >= 0.0 ? "+" : "", @sprintf("%.2f", diff * 100), "%")
end

function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        value, units = t / 1e9, "s"
    end
    return string(@sprintf("%.3f", value), " ", units)
end

function prettymemory(b)
    if b < 1024
        return string(round(Int, b), " bytes")
    elseif b < 1024^2
        value, units = b / 1024, "KiB"
    elseif b < 1024^3
        value, units = b / 1024^2, "MiB"
    else
        value, units = b / 1024^3, "GiB"
    end
    return string(@sprintf("%.2f", value), " ", units)
end

function withtypename(f, io, t)
    needtype = get(io, :typeinfo, Nothing) !== typeof(t)
    if needtype
        print(io, nameof(typeof(t)), '(')
    end
    f()
    if needtype
        print(io, ')')
    end
end

function bindata(sorteddata, nbins, min, max)
    Δ = (max - min) / nbins
    bins = zeros(nbins)
    lastpos = 0
    for i ∈ 1:nbins
        pos = searchsortedlast(sorteddata, min + i * Δ)
        bins[i] = pos - lastpos
        lastpos = pos
    end
    bins
end
bindata(sorteddata, nbins) = bindata(sorteddata, nbins, first(sorteddata), last(sorteddata))

function asciihist(bins, height=1)
    histbars = ['▁', '▂', '▃', '▄', '▅', '▆', '▇', '█']
    if minimum(bins) == 0
        barheights = 2 .+ round.(Int, (height * length(histbars) - 2) * bins ./ maximum(bins))
        barheights[bins .== 0] .= 1
    else
        barheights = 1 .+ round.(Int, (height * length(histbars) - 1) * bins ./ maximum(bins))
    end
    heightmatrix = [min(length(histbars), barheights[b] - (h-1) * length(histbars))
                    for h in height:-1:1, b in 1:length(bins)]
    map(height -> if height < 1; ' ' else histbars[height] end, heightmatrix)
end

_summary(io, t, args...) = withtypename(() -> print(io, args...), io, t)

Base.summary(io::IO, t::Trial) = _summary(io, t, prettytime(time(t)))
Base.summary(io::IO, t::TrialEstimate) = _summary(io, t, prettytime(time(t)))
Base.summary(io::IO, t::TrialRatio) = _summary(io, t, prettypercent(time(t)))
Base.summary(io::IO, t::TrialJudgement) = withtypename(io, t) do
    print(io, prettydiff(time(ratio(t))), " => ")
    printtimejudge(io, t)
end

_show(io, t) =
    if get(io, :compact, true)
        summary(io, t)
    else
        show(io, MIME"text/plain"(), t)
    end

Base.show(io::IO, t::Trial) = _show(io, t)
Base.show(io::IO, t::TrialEstimate) = _show(io, t)
Base.show(io::IO, t::TrialRatio) = _show(io, t)
Base.show(io::IO, t::TrialJudgement) = _show(io, t)

function Base.show(io::IO, m::MIME"text/plain", t::Trial)
    if PREFS.benchmark_output == "fancy"
        show_output(TrialFancyOutputPre(), io, m, t)
    else
        # by default we always use classical version
        show_output(TrialClassicalOutputPre(), io, m, t)
    end

    if PREFS.benchmark_histogram == "fancy"
        show_output(TrialFancyHistogram(), io, m, t)
    else
        # by default we always use classical version
        show_output(TrialClassicalHistogram(), io, m, t)
    end

    if PREFS.benchmark_output == "fancy"
        show_output(TrialFancyOutputPost(), io, m, t)
    else
        # by default we always use classical version
        show_output(TrialClassicalOutputPost(), io, m, t)
    end
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialEstimate)
    println(io, "BenchmarkExt.TrialEstimate: ")
    pad = get(io, :pad, "")
    println(io, pad, "  time:             ", prettytime(time(t)))
    println(io, pad, "  gctime:           ", prettytime(gctime(t)), " (", prettypercent(gctime(t) / time(t)),")")
    println(io, pad, "  memory:           ", prettymemory(memory(t)))
    print(io,   pad, "  allocs:           ", iallocs(t))
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialRatio)
    println(io, "BenchmarkExt.TrialRatio: ")
    pad = get(io, :pad, "")
    println(io, pad, "  time:             ", time(t))
    println(io, pad, "  gctime:           ", gctime(t))
    println(io, pad, "  memory:           ", memory(t))
    print(io,   pad, "  allocs:           ", allocs(t))
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialJudgement)
    println(io, "BenchmarkExt.TrialJudgement: ")
    pad = get(io, :pad, "")
    print(io, pad, "  time:   ", prettydiff(time(ratio(t))), " => ")
    printtimejudge(io, t)
    println(io, " (", prettypercent(params(t).time_tolerance), " tolerance)")
    print(io,   pad, "  memory: ", prettydiff(memory(ratio(t))), " => ")
    printmemoryjudge(io, t)
    println(io, " (", prettypercent(params(t).memory_tolerance), " tolerance)")
end


text:
include("macro.jl")
include("groebner.jl")
include("solve.jl")


text:
module ConfigTests

using ReverseDiff, Test

include(joinpath(dirname(@__FILE__), "../utils.jl"))

issimilar(x::Nothing, y::Nothing) = true
issimilar(x::AbstractArray, y::AbstractArray) = typeof(x) === typeof(y) && size(x) === size(y)
issimilar(x::GradientConfig, y::GradientConfig) = issimilar(x.input, y.input) && x.tape === y.tape
issimilar(x::JacobianConfig, y::JacobianConfig) = issimilar(x.output, y.output) && issimilar(x.input, y.input) && x.tape === y.tape
issimilar(x::Tuple, y::Tuple) = all(map(issimilar, x, y))

############################################################################################

###########
# Config #
###########

tp = InstructionTape()
x, y = rand(3), rand(3, 2)
z = rand(Int, 4)

for Config in (GradientConfig, JacobianConfig)
    cfg = Config(x, tp)
    @test issimilar(cfg.input, track(x, tp))
    @test cfg.tape === tp

    cfg = Config(x, Int, tp)
    @test issimilar(cfg.input, track(x, Int, tp))
    @test cfg.tape === tp

    cfg = Config((x, y), tp)
    @test issimilar(cfg.input, (track(x, tp), track(y, tp)))
    @test cfg.tape === tp

    cfg = Config((x, y), Int, tp)
    @test issimilar(cfg.input, (track(x, Int, tp), track(y, Int, tp)))
    @test cfg.tape === tp
end

cfg = JacobianConfig(z, x, tp)
zt = similar(z, ReverseDiff.TrackedReal{eltype(z),eltype(z),Nothing})
@test issimilar(cfg.input, track(x, eltype(z), tp))
@test issimilar(cfg.output, track!(zt, z, tp))
@test cfg.tape === tp

cfg = JacobianConfig(z, (x, y), tp)
zt = similar(z, ReverseDiff.TrackedReal{eltype(z),eltype(z),Nothing})
@test issimilar(cfg.output, track!(zt, z, tp))
@test issimilar(cfg.input, (track(x, eltype(z), tp), track(y, eltype(z), tp)))
@test cfg.tape === tp

cfg1 = JacobianConfig(z, (x, y), tp)
zt = similar(z, ReverseDiff.TrackedReal{eltype(z),eltype(z),Nothing})
cfg2 = JacobianConfig(DiffResults.JacobianResult(z), (x, y), tp)
@test issimilar(cfg1, cfg2)

##################
# HessianConfig #
##################

gtp, jtp = InstructionTape(), InstructionTape()
x, y = rand(3), rand(Int, 3)

cfg = HessianConfig(x, gtp, jtp)
@test issimilar(cfg.gradient_config, GradientConfig(track(x), gtp))
@test issimilar(cfg.jacobian_config, JacobianConfig(x, jtp))

cfg = HessianConfig(x, Int, gtp, jtp)
@test issimilar(cfg.gradient_config, GradientConfig(track(x, Int), gtp))
@test issimilar(cfg.jacobian_config, JacobianConfig(x, Int, jtp))

cfg = HessianConfig(DiffResults.HessianResult(y), x, gtp, jtp)

# These two act weird: work locally but fail in test mode
#@test issimilar(cfg.gradient_config, GradientConfig(track(x, Int), gtp))
#@test issimilar(cfg.jacobian_config, JacobianConfig(y, x, jtp))

end # module


text:
###
### Constructors
###
###
### Constructor methods for LongSequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

function LongSequence{A}(len::Integer) where {A<:Alphabet}
    if len < 0
        throw(ArgumentError("len must be non-negative"))
    end
    return LongSequence{A}(Vector{UInt64}(undef, seq_data_len(A, len)), 1:convert(Int, len), false)
end

LongSequence(::Type{DNA}) = LongDNASeq()
LongSequence(::Type{RNA}) = LongRNASeq()
LongSequence(::Type{AminoAcid}) = LongAminoAcidSeq()
LongSequence(::Type{Char}) = LongCharSeq()

function LongSequence()
    return LongSequence{VoidAlphabet}(Vector{UInt64}(), 0:-1, false)
end

function LongSequence{A}(s::Union{String, SubString{String}}) where {A<:Alphabet}
    return LongSequence{A}(s, codetype(A()))
end

# Generic method for String/Substring.
function LongSequence{A}(s::Union{String, SubString{String}}, ::AlphabetCode) where {A<:Alphabet}
    len = length(s)
    seq = LongSequence{A}(len)
    return copyto!(seq, 1, s, 1, len)
end

function LongSequence{A}(s::Union{String, SubString{String}}, ::AsciiAlphabet) where {A<:Alphabet}
    v = GC.@preserve s unsafe_wrap(Vector{UInt8}, pointer(s), ncodeunits(s))
    seq = LongSequence{A}(length(v))
    return encode_chunks!(seq, 1, v, 1, length(v))
end

function LongSequence{A}(
        src::Union{AbstractString,AbstractVector},
        startpos::Integer=1,
        stoppos::Integer=length(src)) where {A<:Alphabet}
    len = stoppos - startpos + 1
    seq = LongSequence{A}(len)
    return copyto!(seq, 1, src, startpos, len)
end

# create a subsequence
function LongSequence(other::LongSequence{A}, part::UnitRange{<:Integer}) where {A}
    checkbounds(other, part)
    start = other.part.start + part.start - 1
    stop = start + length(part) - 1
    subseq = LongSequence{A}(other.data, start:stop, true)
    other.shared = true
    return subseq
end

function LongSequence{A}(other::LongSequence{A}, part::UnitRange) where {A}
    return LongSequence(other, part)
end

# Create a 4 bit DNA/RNA sequences from a 2 bit DNA/RNA sequences, and vice-versa.
for (alpha, alphb) in [(DNAAlphabet{4}, DNAAlphabet{2}), # DNA to DNA
                       (DNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{4}, RNAAlphabet{2}), # RNA to RNA
                       (RNAAlphabet{2}, RNAAlphabet{4}),
                       (DNAAlphabet{2}, RNAAlphabet{4}), # DNA to RNA
                       (DNAAlphabet{4}, RNAAlphabet{2}),
                       (RNAAlphabet{4}, DNAAlphabet{2}), # RNA to DNA
                       (RNAAlphabet{2}, DNAAlphabet{4})]

    @eval function (::Type{LongSequence{$alpha}})(seq::LongSequence{$alphb})
        newseq = LongSequence{$alpha}(length(seq))
        for (i, x) in enumerate(seq)
            unsafe_setindex!(newseq, x, i)
        end
        return newseq
    end
end

for (alpha, alphb) in [(DNAAlphabet{2}, RNAAlphabet{2}),
                       (RNAAlphabet{2}, DNAAlphabet{2}),
                       (DNAAlphabet{4}, RNAAlphabet{4}),
                       (RNAAlphabet{4}, DNAAlphabet{4})]

    @eval function (::Type{LongSequence{$alpha}})(seq::LongSequence{$alphb})
        newseq = LongSequence{$alpha}(seq.data, seq.part, true)
        seq.shared = true
        return newseq
    end
end

# Concatenate multiple sequences
function LongSequence{A}(chunks::LongSequence{A}...) where {A}
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = LongSequence{A}(len)
    offset = 1
    for chunk in chunks
        copyto!(seq, offset, chunk, 1, length(chunk))
        offset += length(chunk)
    end
    return seq
end

function Base.repeat(chunk::LongSequence{A}, n::Integer) where {A}
    seq = LongSequence{A}(length(chunk) * n)
    offset = 1
    for _ in 1:n
        copyto!(seq, offset, chunk, 1, length(chunk))
        offset += length(chunk)
    end
    return seq
end

# Concatenation and Base.repeat operators
Base.:*(chunk::LongSequence{A}, chunks::LongSequence{A}...) where {A} =
    LongSequence{A}(chunk, chunks...)
Base.:^(chunk::LongSequence, n::Integer) = repeat(chunk, n)

function Base.similar(seq::LongSequence{A}, len::Integer = length(seq)) where {A}
    return LongSequence{A}(len)
end


text:
@testset "dotproduct" begin
    A = rand(10,5)
    B = rand(20,5)
    d = KernelFunctions.DotProduct()
    @test diag(pairwise(d,A,dims=2)) == [dot(A[:,i],A[:,i]) for i in 1:size(A,2)]
    @test_throws DimensionMismatch d(rand(3),rand(4))
    @test d(3.0,2.0) == 6.0
end


text:
const PyJuliaDictValue_Type__ref = Ref(PyNULL)
PyJuliaDictValue_Type() = begin
    ptr = PyJuliaDictValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.DictValue",
                base = base,
                iter = pyjldict_iter,
                methods = [
                    (name = "keys", flags = Py_METH_NOARGS, meth = pyjldict_keys),
                    (name = "values", flags = Py_METH_NOARGS, meth = pyjldict_values),
                    (name = "items", flags = Py_METH_NOARGS, meth = pyjldict_items),
                    (name = "get", flags = Py_METH_VARARGS, meth = pyjldict_get),
                    (name = "clear", flags = Py_METH_NOARGS, meth = pyjldict_clear),
                    (name = "pop", flags = Py_METH_VARARGS, meth = pyjldict_pop),
                    (name = "popitem", flags = Py_METH_NOARGS, meth = pyjldict_popitem),
                    # (name="update", flags=Py_METH_VARARGS|Py_METH_KEYWORDS, meth=pyjldict_update),
                    (
                        name = "setdefault",
                        flags = Py_METH_VARARGS,
                        meth = pyjldict_setdefault,
                    ),
                ],
                as_sequence = (contains = pyjldict_contains,),
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyMutableMappingABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaDictValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaDictValue_New(x::AbstractDict) = PyJuliaValue_New(PyJuliaDictValue_Type(), x)
PyJuliaValue_From(x::AbstractDict) = PyJuliaDictValue_New(x)

pyjldict_iter(xo::PyPtr) =
    PyJuliaIteratorValue_New(Iterator(keys(PyJuliaValue_GetValue(xo)::AbstractDict)))

pyjldict_keys(xo::PyPtr, ::PyPtr) = PyObject_From(keys(PyJuliaValue_GetValue(xo)))
pyjldict_values(xo::PyPtr, ::PyPtr) = PyObject_From(values(PyJuliaValue_GetValue(xo)))
pyjldict_items(xo::PyPtr, ::PyPtr) =
    PyObject_From(DictPairSet(pairs(PyJuliaValue_GetValue(xo))))

pyjldict_contains(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    r = PyObject_TryConvert(vo, keytype(x))
    r == -1 && return Cint(-1)
    r == 0 && return Cint(0)
    v = takeresult(keytype(x))
    try
        Cint(haskey(x, v))
    catch err
        PyErr_SetJuliaError(err)
        Cint(-1)
    end
end

pyjldict_get(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    ism1(PyArg_CheckNumArgsBetween("get", args, 1, 2)) && return PyNULL
    ko = PyTuple_GetItem(args, 0)
    vo = PyTuple_Size(args) < 2 ? Py_None() : PyTuple_GetItem(args, 1)
    r = PyObject_TryConvert(ko, keytype(x))
    r == -1 && return PyNULL
    r == 0 && (Py_IncRef(vo); return vo)
    k = takeresult(keytype(x))
    try
        if haskey(x, k)
            PyObject_From(x[k])
        else
            Py_IncRef(vo)
            vo
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjldict_setdefault(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    ism1(PyArg_CheckNumArgsBetween("get", args, 1, 2)) && return PyNULL
    ko = PyTuple_GetItem(args, 0)
    vo = PyTuple_Size(args) < 2 ? Py_None() : PyTuple_GetItem(args, 1)
    r = PyObject_TryConvert(ko, keytype(x))
    r == -1 && return PyNULL
    r == 0 && (Py_IncRef(vo); return vo)
    k = takeresult(keytype(x))
    try
        if haskey(x, k)
            PyObject_From(x[k])
        else
            r = PyObject_Convert(vo, valtype(x))
            r == -1 && return PyNULL
            x[k] = takeresult(valtype(x))
            Py_IncRef(vo)
            vo
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjldict_clear(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    try
        empty!(x)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjldict_popitem(xo::PyPtr, ::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    try
        if isempty(x)
            PyErr_SetString(PyExc_KeyError(), "pop empty dictionary")
            PyNULL
        else
            k, v = pop!(x)
            PyTuple_From((k, v))
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjldict_pop(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractDict
    ism1(PyArg_CheckNumArgsBetween("pop", args, 1, 2)) && return PyNULL
    ko = PyTuple_GetItem(args, 0)
    vo = PyTuple_Size(args) == 2 ? PyTuple_GetItem(args, 1) : PyNULL
    r = PyObject_TryConvert(PyTuple_GetItem(args, 0), keytype(x))
    r == -1 && return PyNULL
    r == 0 &&
        (isnull(vo) ? (PyErr_SetObject(PyExc_KeyError(), ko)) : (Py_IncRef(vo); vo))
    k = takeresult(keytype(x))
    try
        if haskey(x, k)
            PyObject_From(pop!(x, k))
        elseif isnull(vo)
            PyErr_SetObject(PyExc_KeyError(), ko)
            PyNULL
        else
            Py_IncRef(vo)
            vo
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

struct DictPairSet{K,V,T<:AbstractDict{K,V}} <: AbstractSet{Tuple{K,V}}
    dict::T
end
Base.length(x::DictPairSet) = length(x.dict)
Base.iterate(x::DictPairSet) =
    (r = iterate(x.dict); r === nothing ? nothing : (Tuple(r[1]), r[2]))
Base.iterate(x::DictPairSet, st) =
    (r = iterate(x.dict, st); r === nothing ? nothing : (Tuple(r[1]), r[2]))
Base.in(v, x::DictPairSet) = v in x.dict


text:
using TopOptDash
using Test

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


text:
using Base.Test, DiffEqBase, SimpleTraits

f123(x,y) = 1
f123(::Val{:jac}, x) = 2

@test istrait(HasJac{typeof(f123)})
@code_llvm has_jac(f123)

g123(x,y) = 1
g123{T}(::Val{:jac}, x::T) = 2
@traitimpl HasJac{typeof(g123)}
@code_llvm has_jac(g123)

h123(x,y) = 1
h123{T}(x::T) = 2

immutable T123 end
(::T123)(a) = 1
(::T123)(::Val{:jac}, a) = 1
t123 = T123()
@code_llvm has_jac(t123)

immutable G123 end
(::G123)(a) = 1
(::G123)(b, a) = 1

@test has_jac(f123)
@test has_jac(g123)
@test !has_jac(h123)
@test has_jac(T123())
@test !has_jac(G123())

@traitfn tfn(f::::HasJac) = true
@traitfn tfn(f::::(!HasJac)) = false

@traitfn testing(f::::HasJac)   =2
@traitfn testing(f::::!(HasJac))=3.0

@inferred testing(f123)
@inferred testing(g123)
@inferred testing(h123)

@test tfn(f123)
@test tfn(g123)
@test !tfn(h123)
@test tfn(T123())
@test !tfn(G123())
@test !tfn(T123)  # NOTE this is inconsistent as has_jac(T123)==true
@test !tfn(G123)

using ParameterizedFunctions

type  LotkaVolterra <: ParameterizedFunction
          a::Float64
          b::Float64
end
(p::LotkaVolterra)(t,u,du) = begin
        du[1] = p.a * u[1] - p.b * u[1]*u[2]
        du[2] = -3 * u[2] + u[1]*u[2]
end
(p::LotkaVolterra)(::Val{:jac}, t, u, J) = 1

lv = LotkaVolterra(0.0,0.0)

@test has_jac(lv)
@test tfn(lv)


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "image1"
    help = "The first comparison image"
    "config1"
    help = "a YAML configuration file for the first image"
    "image2"
    help = "The second comparison image"
    "config2"
    help = "a YAML configuration file for the second image"
    "--per"
    help = "Percentage to scale the residual."
    default = 0.1 # 10% is max
    "--log"
    help = "Should the intensity be plotted in log scale?"
    action = :store_true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config1"]))
config2 = YAML.load(open(parsed_args["config2"]))

# Load two separate high resolution images and compare them, computing the difference and the difference in lnprob.

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.gridding

import PyPlot.plt
using LaTeXStrings

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

pars = convert_dict(config["parameters"], config["model"])
pars2 = convert_dict(config2["parameters"], config2["model"])

img1 = imread(parsed_args["image1"])
img2 = imread(parsed_args["image2"])

nchan = length(img1.lams)
# Make sure we have the same number of channels to compare images.
@assert nchan == length(img2.lams) "Images have a different number of channels!"

# Convert raw images to the appropriate distance
skimg1 = imToSky(img1, pars.dpc)
skimg2 = imToSky(img2, pars2.dpc)
skresid = skimg1 - skimg2

# Store maximums for plotting
max_img = maxabs(skimg1.data)
if maxabs(skimg2.data) > max_img
    max_img = maxabs(skimg2.data)
end
max_resid = maxabs(skresid.data)

# Apply the gridding correction function before doing the FFT
# No shift needed, since we will shift the resampled visibilities
# The skimg1, skimg2 are for plotting.
skimg01 = corrfun(skimg1)
skimg02 = corrfun(skimg2)

# To store the FFT for each channel
vis_fft1 = Array(FullModelVis, nchan)
vis_fft2 = Array(FullModelVis, nchan)
vis_resid = Array(FullModelVis, nchan)

# To store the computed lnprob for each channel
lnprobs1 = Array(Float64, nchan)
lnprobs2 = Array(Float64, nchan)
lnprobs_diff = Array(Float64, nchan)

# Load the data file so we can use it to compute lnprob
dvarr = DataVis(config["data_file"])
for dset in dvarr
    # Conjugation is necessary for the SMA and ALMA
    visibilities.conj!(dset) # Swap UV convention
end

half_pix = config["size_arcsec"] / (2 * config["npix"])

# Fourier transform each channel, while also noting the max scale values
for i=1:nchan
    dv = dvarr[i]
    # FFT the appropriate image channel
    vis_fft1[i] = transform(skimg01, i)
    vis_fft2[i] = transform(skimg02, i)
    vis_resid[i] = vis_fft1[i] - vis_fft2[i]

    mvis1 = ModelVis(dv, vis_fft1[i])
    mvis2 = ModelVis(dv, vis_fft2[i])

    # Apply the phase shift here
    phase_shift!(mvis1, pars.mu_RA + half_pix, pars.mu_DEC - half_pix)
    phase_shift!(mvis2, pars2.mu_RA + half_pix, pars2.mu_DEC - half_pix)

    # Calculate the likelihood between these two
    lnprobs1[i] = lnprob(dv, mvis1)
    lnprobs2[i] = lnprob(dv, mvis2)
    lnprobs_diff[i] = lnprobs1[i] - lnprobs2[i]
end

# Store the max and min values for each type of plot
max_real = maxabs(real(vis_fft1[1].VV))
max_real_resid = maxabs(real(vis_resid[1].VV))
max_imag = maxabs(imag(vis_fft1[1].VV))
max_imag_resid = maxabs(imag(vis_resid[1].VV))

for i=1:nchan

    max_real_temp = maxabs(real(vis_fft1[i].VV))
    if max_real_temp > max_real
        max_real = max_real_temp
    end

    max_real_temp = maxabs(real(vis_fft2[i].VV))
    if max_real_temp > max_real
        max_real = max_real_temp
    end

    max_real_temp = maxabs(real(vis_resid[i].VV))
    if max_real_temp > max_real_resid
        max_real_resid = max_real_temp
    end

    max_imag_temp = maxabs(imag(vis_fft1[i].VV))
    if max_imag_temp > max_imag
        max_imag = max_imag_temp
    end

    max_imag_temp = maxabs(imag(vis_fft2[i].VV))
    if max_imag_temp > max_imag
        max_imag = max_imag_temp
    end

    max_imag_temp = maxabs(imag(vis_resid[i].VV))
    if max_imag_temp > max_imag_resid
        max_imag_resid = max_imag_temp
    end

end

# Create the normalize instances for each
norm_img = PyPlot.matplotlib[:colors][:Normalize](0, max_img)
norm_resid = PyPlot.matplotlib[:colors][:Normalize](-max_resid, max_resid)

norm_real = PyPlot.matplotlib[:colors][:Normalize](-max_real, max_real)
norm_real_resid = PyPlot.matplotlib[:colors][:Normalize](-max_real_resid, max_real_resid)

norm_imag = PyPlot.matplotlib[:colors][:Normalize](-max_imag, max_imag)
norm_imag_resid = PyPlot.matplotlib[:colors][:Normalize](-max_imag_resid, max_imag_resid)

# Now go through and plot everything
fig, ax = plt[:subplots](nrows=10, ncols=nchan, figsize=(2.0 * nchan, 24.))

# Image needs to be flipped along RA dimension
ext_img = (skimg1.ra[end], skimg1.ra[1], skimg1.dec[1], skimg1.dec[end])

# extent in Fourier domain
ext_vis = (vis_fft1[1].uu[1], vis_fft1[1].uu[end], vis_fft1[1].vv[1], vis_fft1[1].vv[end])

cm = plt[:get_cmap]("PuBu")
# Reversed color map
cm_resid = plt[:get_cmap]("bwr_r")

for i=1:nchan
    # Images
    ax[1, i][:imshow](flipdim(skimg1.data[:,:,i], 2), interpolation="none", origin="lower", cmap=cm, extent=ext_img, norm=norm_img)
    ax[2, i][:imshow](flipdim(skimg2.data[:,:,i], 2), interpolation="none", origin="lower", cmap=cm, extent=ext_img, norm=norm_img)
    im = ax[3, i][:imshow](flipdim(skresid.data[:,:,i], 2), interpolation="none", origin="lower", cmap=cm_resid, extent=ext_img, norm=norm_resid)

    if i==nchan
        cax = fig[:add_axes]([0.95, 0.7, 0.01, 0.07])
        cbar = fig[:colorbar](im, cax=cax)
        # cbar.ax.tick_params(labelsize=6)
    end

    ## Fourier Domain

    # Real
    ax[4, i][:imshow](real(vis_fft1[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_real)
    ax[5, i][:imshow](real(vis_fft2[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_real)
    im = ax[6, i][:imshow](real(vis_resid[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_real_resid)

    if i==nchan
        cax = fig[:add_axes]([0.95, 0.43, 0.01, 0.07])
        cbar = fig[:colorbar](im, cax=cax)
        # cbar.ax.tick_params(labelsize=6)
    end

    # Imag
    ax[7, i][:imshow](imag(vis_fft1[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_imag)
    ax[8, i][:imshow](imag(vis_fft2[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_imag)
    im = ax[9, i][:imshow](imag(vis_resid[i].VV), extent=ext_vis, interpolation="none", origin="lower", cmap=cm_resid, norm=norm_imag_resid)

    if i==nchan
        cax = fig[:add_axes]([0.95, 0.15, 0.01, 0.07])
        cbar = fig[:colorbar](im, cax=cax)
        # cbar.ax.tick_params(labelsize=6)
    end

    ax[9, i][:annotate](@sprintf("%.2f", lnprobs_diff[i]), (0.1, 0.2), xycoords="axes fraction", size=8)

    dv = dvarr[i]
    ax[10,i][:plot](dv.uu, dv.vv, "k.", ms=0.2)
    ax[10,i][:set_xlim](vis_fft1[1].uu[1], vis_fft1[1].uu[end])
    ax[10,i][:set_ylim](vis_fft1[1].vv[1], vis_fft1[1].vv[end])

    # Remove tick labels from all subsequent columns.
    if i==1
        # Do nothing
    else
        for j=1:10
            ax[j, i][:xaxis][:set_ticklabels]([])
            ax[j, i][:yaxis][:set_ticklabels]([])
        end
    end

    # for j=1:10
    #     if (j==1 && i==1) || (j==4 && i==1) || (j==10 && i==1)
    #         # Do nothing
    #     else
    #         ax[j, i][:xaxis][:set_ticklabels]([])
    #         ax[j, i][:yaxis][:set_ticklabels]([])
    #     end
    # end
end

for j=1:3
    ax[j, 1][:set_xlabel](L"$\Delta \alpha$ ['']")
    ax[j, 1][:set_ylabel](L"$\Delta \delta$ ['']")
end

for j=4:10
    ax[j, 1][:set_ylabel](L"UU $[k\lambda]$")
    ax[j, 1][:set_xlabel](L"VV $[k\lambda]$")
end

name1 = config["name"]
name2 = config2["name"]

ax[1, 1][:annotate]("image $name1", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[2, 1][:annotate]("image $name2", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[3, 1][:annotate]("image $name1 - $name2", (0.1, 0.8), xycoords="axes fraction", size=8)

ax[4, 1][:annotate]("real $name1", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[5, 1][:annotate]("real $name2", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[6, 1][:annotate]("real $name1 - $name2", (0.1, 0.8), xycoords="axes fraction", size=8)

ax[7, 1][:annotate]("imag $name1", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[8, 1][:annotate]("imag $name2", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[9, 1][:annotate]("imag $name1 - $name2", (0.1, 0.8), xycoords="axes fraction", size=8)
ax[9, 1][:annotate]("lnprob diff", (0.1, 0.4), xycoords="axes fraction", size=8)
ax[10, 1][:annotate]("baselines", (0.1, 0.8), xycoords="axes fraction", size=8)

fig[:subplots_adjust](hspace=0.3, wspace=0.001, top=0.95, bottom=0.05, left=0.1, right=0.9)
fig[:savefig]("chmaps_compare.png", dpi=300)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using DynamicPPL
using DynamicPPL: apply_dotted, getargs_dottilde, getargs_tilde

using Test

@testset "apply_dotted" begin
    # Some things that are not expressions.
    @test apply_dotted(:x) === :x
    @test apply_dotted(1.0) === 1.0
    @test apply_dotted([1.0, 2.0, 4.0]) == [1.0, 2.0, 4.0]

    # Some expressions.
    @test apply_dotted(:(x ~ Normal(μ, σ))) == :(x ~ Normal(μ, σ))
    @test apply_dotted(:((.~)(x, Normal(μ, σ)))) == :((.~)(x, Normal(μ, σ)))
    @test apply_dotted(:((~).(x, Normal(μ, σ)))) == :((~).(x, Normal(μ, σ)))
    @test apply_dotted(:(@. x ~ Normal(μ, σ))) == :((~).(x, Normal.(μ, σ)))
    @test apply_dotted(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) ==
        :((~).(x, Normal.(μ, sqrt(v))))
    @test apply_dotted(:(@~ Normal.(μ, σ))) == :(@~ Normal.(μ, σ))
end

@testset "getargs_dottilde" begin
    # Some things that are not expressions.
    @test getargs_dottilde(:x) === nothing
    @test getargs_dottilde(1.0) === nothing
    @test getargs_dottilde([1.0, 2.0, 4.0]) === nothing

    # Some expressions.
    @test getargs_dottilde(:(x ~ Normal(μ, σ))) === nothing
    @test getargs_dottilde(:((.~)(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
    @test getargs_dottilde(:((~).(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
    @test getargs_dottilde(:(@. x ~ Normal(μ, σ))) === nothing
    @test getargs_dottilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
    @test getargs_dottilde(:(@~ Normal.(μ, σ))) === nothing
end

@testset "getargs_tilde" begin
    # Some things that are not expressions.
    @test getargs_tilde(:x) === nothing
    @test getargs_tilde(1.0) === nothing
    @test getargs_tilde([1.0, 2.0, 4.0]) === nothing

    # Some expressions.
    @test getargs_tilde(:(x ~ Normal(μ, σ))) == (:x, :(Normal(μ, σ)))
    @test getargs_tilde(:((.~)(x, Normal(μ, σ)))) === nothing
    @test getargs_tilde(:((~).(x, Normal(μ, σ)))) === nothing
    @test getargs_tilde(:(@. x ~ Normal(μ, σ))) === nothing
    @test getargs_tilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
    @test getargs_tilde(:(@~ Normal.(μ, σ))) === nothing
end


text:
function safe_symbol(cur_reactor::AbstractReactor, cur_field::Symbol)
  cur_value = getfield(cur_reactor, cur_field)

  isa(cur_value, SymEngine.Basic) || return cur_value

  isempty(free_symbols(cur_value)) && return N(cur_value)

  return symbols(cur_field)
end