text:
using Base.Ryu

const maxMantissa = (UInt64(1) << 53) - 1
todouble(sign, exp, mant) = Core.bitcast(Float64, (UInt64(sign) << 63) | (UInt64(exp) << 52) | (UInt64(mant)))

@testset "Ryu" begin

@testset "Float64" begin

@testset "Basic" begin
    @test Ryu.writeshortest(0.0) == "0.0"
    @test Ryu.writeshortest(-0.0) == "-0.0"
    @test Ryu.writeshortest(1.0) == "1.0"
    @test Ryu.writeshortest(-1.0) == "-1.0"
    @test Ryu.writeshortest(NaN) == "NaN"
    @test Ryu.writeshortest(Inf) == "Inf"
    @test Ryu.writeshortest(-Inf) == "-Inf"
end

@testset "SwitchToSubnormal" begin
    @test "2.2250738585072014e-308" == Ryu.writeshortest(2.2250738585072014e-308)
end

@testset "MinAndMax" begin
    @test "1.7976931348623157e308" == Ryu.writeshortest(Core.bitcast(Float64, 0x7fefffffffffffff))
    @test "5.0e-324" == Ryu.writeshortest(Core.bitcast(Float64, Int64(1)))
end

@testset "LotsOfTrailingZeros" begin
    @test "2.9802322387695312e-8" == Ryu.writeshortest(2.98023223876953125e-8)
end

@testset "Regression" begin
    @test "-2.109808898695963e16" == Ryu.writeshortest(-2.109808898695963e16)
    @test "4.940656e-318" == Ryu.writeshortest(4.940656e-318)
    @test "1.18575755e-316" == Ryu.writeshortest(1.18575755e-316)
    @test "2.989102097996e-312" == Ryu.writeshortest(2.989102097996e-312)
    @test "9.0608011534336e15" == Ryu.writeshortest(9.0608011534336e15)
    @test "4.708356024711512e18" == Ryu.writeshortest(4.708356024711512e18)
    @test "9.409340012568248e18" == Ryu.writeshortest(9.409340012568248e18)
    @test "1.2345678" == Ryu.writeshortest(1.2345678)
end

@testset "LooksLikePow5" begin
    # These numbers have a mantissa that is a multiple of the largest power of 5 that fits,
    # and an exponent that causes the computation for q to result in 22, which is a corner
    # case for Ryu.
    @test "5.764607523034235e39" == Ryu.writeshortest(Core.bitcast(Float64, 0x4830F0CF064DD592))
    @test "1.152921504606847e40" == Ryu.writeshortest(Core.bitcast(Float64, 0x4840F0CF064DD592))
    @test "2.305843009213694e40" == Ryu.writeshortest(Core.bitcast(Float64, 0x4850F0CF064DD592))
end

@testset "OutputLength" begin
    @test "1.0" == Ryu.writeshortest(1.0) # already tested in Basic
    @test "1.2" == Ryu.writeshortest(1.2)
    @test "1.23" == Ryu.writeshortest(1.23)
    @test "1.234" == Ryu.writeshortest(1.234)
    @test "1.2345" == Ryu.writeshortest(1.2345)
    @test "1.23456" == Ryu.writeshortest(1.23456)
    @test "1.234567" == Ryu.writeshortest(1.234567)
    @test "1.2345678" == Ryu.writeshortest(1.2345678) # already tested in Regressi
    @test "1.23456789" == Ryu.writeshortest(1.23456789)
    @test "1.234567895" == Ryu.writeshortest(1.234567895) # 1.234567890 would be trimm
    @test "1.2345678901" == Ryu.writeshortest(1.2345678901)
    @test "1.23456789012" == Ryu.writeshortest(1.23456789012)
    @test "1.234567890123" == Ryu.writeshortest(1.234567890123)
    @test "1.2345678901234" == Ryu.writeshortest(1.2345678901234)
    @test "1.23456789012345" == Ryu.writeshortest(1.23456789012345)
    @test "1.234567890123456" == Ryu.writeshortest(1.234567890123456)
    @test "1.2345678901234567" == Ryu.writeshortest(1.2345678901234567)

  # Test 32-bit chunking
    @test "4.294967294" == Ryu.writeshortest(4.294967294) # 2^32 -
    @test "4.294967295" == Ryu.writeshortest(4.294967295) # 2^32 -
    @test "4.294967296" == Ryu.writeshortest(4.294967296) # 2^
    @test "4.294967297" == Ryu.writeshortest(4.294967297) # 2^32 +
    @test "4.294967298" == Ryu.writeshortest(4.294967298) # 2^32 +
end

# Test min, max shift values in shiftright128
@testset "MinMaxShift" begin
    # 32-bit opt-size=0:  49 <= dist <= 50
    # 32-bit opt-size=1:  30 <= dist <= 50
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  30 <= dist <= 50
    @test "1.7800590868057611e-307" == Ryu.writeshortest(todouble(false, 4, 0))
    # 32-bit opt-size=0:  49 <= dist <= 49
    # 32-bit opt-size=1:  28 <= dist <= 49
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  28 <= dist <= 50
    @test "2.8480945388892175e-306" == Ryu.writeshortest(todouble(false, 6, maxMantissa))
    # 32-bit opt-size=0:  52 <= dist <= 53
    # 32-bit opt-size=1:   2 <= dist <= 53
    # 64-bit opt-size=0:  53 <= dist <= 53
    # 64-bit opt-size=1:   2 <= dist <= 53
    @test "2.446494580089078e-296" == Ryu.writeshortest(todouble(false, 41, 0))
    # 32-bit opt-size=0:  52 <= dist <= 52
    # 32-bit opt-size=1:   2 <= dist <= 52
    # 64-bit opt-size=0:  53 <= dist <= 53
    # 64-bit opt-size=1:   2 <= dist <= 53
    @test "4.8929891601781557e-296" == Ryu.writeshortest(todouble(false, 40, maxMantissa))

    # 32-bit opt-size=0:  57 <= dist <= 58
    # 32-bit opt-size=1:  57 <= dist <= 58
    # 64-bit opt-size=0:  58 <= dist <= 58
    # 64-bit opt-size=1:  58 <= dist <= 58
    @test "1.8014398509481984e16" == Ryu.writeshortest(todouble(false, 1077, 0))
    # 32-bit opt-size=0:  57 <= dist <= 57
    # 32-bit opt-size=1:  57 <= dist <= 57
    # 64-bit opt-size=0:  58 <= dist <= 58
    # 64-bit opt-size=1:  58 <= dist <= 58
    @test "3.6028797018963964e16" == Ryu.writeshortest(todouble(false, 1076, maxMantissa))
    # 32-bit opt-size=0:  51 <= dist <= 52
    # 32-bit opt-size=1:  51 <= dist <= 59
    # 64-bit opt-size=0:  52 <= dist <= 52
    # 64-bit opt-size=1:  52 <= dist <= 59
    @test "2.900835519859558e-216" == Ryu.writeshortest(todouble(false, 307, 0))
    # 32-bit opt-size=0:  51 <= dist <= 51
    # 32-bit opt-size=1:  51 <= dist <= 59
    # 64-bit opt-size=0:  52 <= dist <= 52
    # 64-bit opt-size=1:  52 <= dist <= 59
    @test "5.801671039719115e-216" == Ryu.writeshortest(todouble(false, 306, maxMantissa))

    # https:#github.com/ulfjack/ryu/commit/19e44d16d80236f5de25800f56d82606d1be00b9#commitcomment-30146483
    # 32-bit opt-size=0:  49 <= dist <= 49
    # 32-bit opt-size=1:  44 <= dist <= 49
    # 64-bit opt-size=0:  50 <= dist <= 50
    # 64-bit opt-size=1:  44 <= dist <= 50
    @test "3.196104012172126e-27" == Ryu.writeshortest(todouble(false, 934, 0x000FA7161A4D6E0C))
end

@testset "SmallIntegers" begin
    @test "9.007199254740991e15" == Ryu.writeshortest(9007199254740991.0)
    @test "9.007199254740992e15" == Ryu.writeshortest(9007199254740992.0)

    @test "1.0" == Ryu.writeshortest(1.0e+0)
    @test "12.0" == Ryu.writeshortest(1.2e+1)
    @test "123.0" == Ryu.writeshortest(1.23e+2)
    @test "1234.0" == Ryu.writeshortest(1.234e+3)
    @test "12345.0" == Ryu.writeshortest(1.2345e+4)
    @test "123456.0" == Ryu.writeshortest(1.23456e+5)
    @test "1.234567e6" == Ryu.writeshortest(1.234567e+6)
    @test "1.2345678e7" == Ryu.writeshortest(1.2345678e+7)
    @test "1.23456789e8" == Ryu.writeshortest(1.23456789e+8)
    @test "1.23456789e9" == Ryu.writeshortest(1.23456789e+9)
    @test "1.234567895e9" == Ryu.writeshortest(1.234567895e+9)
    @test "1.2345678901e10" == Ryu.writeshortest(1.2345678901e+10)
    @test "1.23456789012e11" == Ryu.writeshortest(1.23456789012e+11)
    @test "1.234567890123e12" == Ryu.writeshortest(1.234567890123e+12)
    @test "1.2345678901234e13" == Ryu.writeshortest(1.2345678901234e+13)
    @test "1.23456789012345e14" == Ryu.writeshortest(1.23456789012345e+14)
    @test "1.234567890123456e15" == Ryu.writeshortest(1.234567890123456e+15)

  # 10^i
    @test "1.0" == Ryu.writeshortest(1.0e+0)
    @test "10.0" == Ryu.writeshortest(1.0e+1)
    @test "100.0" == Ryu.writeshortest(1.0e+2)
    @test "1000.0" == Ryu.writeshortest(1.0e+3)
    @test "10000.0" == Ryu.writeshortest(1.0e+4)
    @test "100000.0" == Ryu.writeshortest(1.0e+5)
    @test "1.0e6" == Ryu.writeshortest(1.0e+6)
    @test "1.0e7" == Ryu.writeshortest(1.0e+7)
    @test "1.0e8" == Ryu.writeshortest(1.0e+8)
    @test "1.0e9" == Ryu.writeshortest(1.0e+9)
    @test "1.0e10" == Ryu.writeshortest(1.0e+10)
    @test "1.0e11" == Ryu.writeshortest(1.0e+11)
    @test "1.0e12" == Ryu.writeshortest(1.0e+12)
    @test "1.0e13" == Ryu.writeshortest(1.0e+13)
    @test "1.0e14" == Ryu.writeshortest(1.0e+14)
    @test "1.0e15" == Ryu.writeshortest(1.0e+15)

  # 10^15 + 10^i
    @test "1.000000000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+0)
    @test "1.00000000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+1)
    @test "1.0000000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+2)
    @test "1.000000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+3)
    @test "1.00000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+4)
    @test "1.0000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+5)
    @test "1.000000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+6)
    @test "1.00000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+7)
    @test "1.0000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+8)
    @test "1.000001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+9)
    @test "1.00001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+10)
    @test "1.0001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+11)
    @test "1.001e15" == Ryu.writeshortest(1.0e+15 + 1.0e+12)
    @test "1.01e15" == Ryu.writeshortest(1.0e+15 + 1.0e+13)
    @test "1.1e15" == Ryu.writeshortest(1.0e+15 + 1.0e+14)

  # Largest power of 2 <= 10^(i+1)
    @test "8.0" == Ryu.writeshortest(8.0)
    @test "64.0" == Ryu.writeshortest(64.0)
    @test "512.0" == Ryu.writeshortest(512.0)
    @test "8192.0" == Ryu.writeshortest(8192.0)
    @test "65536.0" == Ryu.writeshortest(65536.0)
    @test "524288.0" == Ryu.writeshortest(524288.0)
    @test "8.388608e6" == Ryu.writeshortest(8388608.0)
    @test "6.7108864e7" == Ryu.writeshortest(67108864.0)
    @test "5.36870912e8" == Ryu.writeshortest(536870912.0)
    @test "8.589934592e9" == Ryu.writeshortest(8589934592.0)
    @test "6.8719476736e10" == Ryu.writeshortest(68719476736.0)
    @test "5.49755813888e11" == Ryu.writeshortest(549755813888.0)
    @test "8.796093022208e12" == Ryu.writeshortest(8796093022208.0)
    @test "7.0368744177664e13" == Ryu.writeshortest(70368744177664.0)
    @test "5.62949953421312e14" == Ryu.writeshortest(562949953421312.0)
    @test "9.007199254740992e15" == Ryu.writeshortest(9007199254740992.0)

  # 1000 * (Largest power of 2 <= 10^(i+1))
    @test "8000.0" == Ryu.writeshortest(8.0e+3)
    @test "64000.0" == Ryu.writeshortest(64.0e+3)
    @test "512000.0" == Ryu.writeshortest(512.0e+3)
    @test "8.192e6" == Ryu.writeshortest(8192.0e+3)
    @test "6.5536e7" == Ryu.writeshortest(65536.0e+3)
    @test "5.24288e8" == Ryu.writeshortest(524288.0e+3)
    @test "8.388608e9" == Ryu.writeshortest(8388608.0e+3)
    @test "6.7108864e10" == Ryu.writeshortest(67108864.0e+3)
    @test "5.36870912e11" == Ryu.writeshortest(536870912.0e+3)
    @test "8.589934592e12" == Ryu.writeshortest(8589934592.0e+3)
    @test "6.8719476736e13" == Ryu.writeshortest(68719476736.0e+3)
    @test "5.49755813888e14" == Ryu.writeshortest(549755813888.0e+3)
    @test "8.796093022208e15" == Ryu.writeshortest(8796093022208.0e+3)
end

end # Float64

@testset "Float32" begin

@testset "Basic" begin
    @test "0.0" == Ryu.writeshortest(Float32(0.0))
    @test "-0.0" == Ryu.writeshortest(Float32(-0.0))
    @test "1.0" == Ryu.writeshortest(Float32(1.0))
    @test "-1.0" == Ryu.writeshortest(Float32(-1.0))
    @test "NaN" == Ryu.writeshortest(Float32(NaN))
    @test "Inf" == Ryu.writeshortest(Float32(Inf))
    @test "-Inf" == Ryu.writeshortest(Float32(-Inf))
end

@testset "SwitchToSubnormal" begin
    @test "1.1754944e-38" == Ryu.writeshortest(1.1754944f-38)
end

@testset "MinAndMax" begin
    @test "3.4028235e38" == Ryu.writeshortest(Core.bitcast(Float32, 0x7f7fffff))
    @test "1.0e-45" == Ryu.writeshortest(Core.bitcast(Float32, Int32(1)))
end

# Check that we return the exact boundary if it is the shortest
# representation, but only if the original floating point number is even.
@testset "BoundaryRoundeven" begin
    @test "3.355445e7" == Ryu.writeshortest(3.355445f7)
    @test "9.0e9" == Ryu.writeshortest(8.999999f9)
    @test "3.436672e10" == Ryu.writeshortest(3.4366717f10)
end

# If the exact value is exactly halfway between two shortest representations,
# then we round to even. It seems like this only makes a difference if the
# last two digits are ...2|5 or ...7|5, and we cut off the 5.
@testset "exactValueRoundeven" begin
    @test "305404.12" == Ryu.writeshortest(3.0540412f5)
    @test "8099.0312" == Ryu.writeshortest(8.0990312f3)
end

@testset "LotsOfTrailingZeros" begin
    # Pattern for the first test: 00111001100000000000000000000000
    @test "0.00024414062" == Ryu.writeshortest(2.4414062f-4)
    @test "0.0024414062" == Ryu.writeshortest(2.4414062f-3)
    @test "0.0043945312" == Ryu.writeshortest(4.3945312f-3)
    @test "0.0063476562" == Ryu.writeshortest(6.3476562f-3)
end

@testset "Regression" begin
    @test "4.7223665e21" == Ryu.writeshortest(4.7223665f21)
    @test "8.388608e6" == Ryu.writeshortest(8388608f0)
    @test "1.6777216e7" == Ryu.writeshortest(1.6777216f7)
    @test "3.3554436e7" == Ryu.writeshortest(3.3554436f7)
    @test "6.7131496e7" == Ryu.writeshortest(6.7131496f7)
    @test "1.9310392e-38" == Ryu.writeshortest(1.9310392f-38)
    @test "-2.47e-43" == Ryu.writeshortest(-2.47f-43)
    @test "1.993244e-38" == Ryu.writeshortest(1.993244f-38)
    @test "4103.9004" == Ryu.writeshortest(4103.9003f0)
    @test "5.3399997e9" == Ryu.writeshortest(5.3399997f9)
    @test "6.0898e-39" == Ryu.writeshortest(6.0898f-39)
    @test "0.0010310042" == Ryu.writeshortest(0.0010310042f0)
    @test "2.882326e17" == Ryu.writeshortest(2.8823261f17)
    @test "7.038531e-26" == Ryu.writeshortest(7.0385309f-26)
    @test "9.223404e17" == Ryu.writeshortest(9.2234038f17)
    @test "6.710887e7" == Ryu.writeshortest(6.7108872f7)
    @test "1.0e-44" == Ryu.writeshortest(1.0f-44)
    @test "2.816025e14" == Ryu.writeshortest(2.816025f14)
    @test "9.223372e18" == Ryu.writeshortest(9.223372f18)
    @test "1.5846086e29" == Ryu.writeshortest(1.5846085f29)
    @test "1.1811161e19" == Ryu.writeshortest(1.1811161f19)
    @test "5.368709e18" == Ryu.writeshortest(5.368709f18)
    @test "4.6143166e18" == Ryu.writeshortest(4.6143165f18)
    @test "0.007812537" == Ryu.writeshortest(0.007812537f0)
    @test "1.0e-45" == Ryu.writeshortest(1.4f-45)
    @test "1.18697725e20" == Ryu.writeshortest(1.18697724f20)
    @test "1.00014165e-36" == Ryu.writeshortest(1.00014165f-36)
    @test "200.0" == Ryu.writeshortest(200f0)
    @test "3.3554432e7" == Ryu.writeshortest(3.3554432f7)
end

@testset "LooksLikePow5" begin
    # These numbers have a mantissa that is the largest power of 5 that fits,
    # and an exponent that causes the computation for q to result in 10, which is a corner
    # case for Ryu.
    @test "6.7108864e17" == Ryu.writeshortest(Core.bitcast(Float32, 0x5D1502F9))
    @test "1.3421773e18" == Ryu.writeshortest(Core.bitcast(Float32, 0x5D9502F9))
    @test "2.6843546e18" == Ryu.writeshortest(Core.bitcast(Float32, 0x5E1502F9))
end

@testset "OutputLength" begin
    @test "1.0" == Ryu.writeshortest(Float32(1.0))
    @test "1.2" == Ryu.writeshortest(Float32(1.2))
    @test "1.23" == Ryu.writeshortest(Float32(1.23))
    @test "1.234" == Ryu.writeshortest(Float32(1.234))
    @test "1.2345" == Ryu.writeshortest(Float32(1.2345))
    @test "1.23456" == Ryu.writeshortest(Float32(1.23456))
    @test "1.234567" == Ryu.writeshortest(Float32(1.234567))
    @test "1.2345678" == Ryu.writeshortest(Float32(1.2345678))
    @test "1.23456735e-36" == Ryu.writeshortest(Float32(1.23456735e-36))
end

end # Float32

@testset "Float16" begin

@testset "Basic" begin
    @test "0.0" == Ryu.writeshortest(Float16(0.0))
    @test "-0.0" == Ryu.writeshortest(Float16(-0.0))
    @test "1.0" == Ryu.writeshortest(Float16(1.0))
    @test "-1.0" == Ryu.writeshortest(Float16(-1.0))
    @test "NaN" == Ryu.writeshortest(Float16(NaN))
    @test "Inf" == Ryu.writeshortest(Float16(Inf))
    @test "-Inf" == Ryu.writeshortest(Float16(-Inf))
end

let x=floatmin(Float16)
    while x <= floatmax(Float16)
        @test parse(Float16, Ryu.writeshortest(x)) == x
        x = nextfloat(x)
    end
end

# function testfloats(T)
#     x = floatmin(T)
#     i = 0
#     fails = 0
#     success = 0
#     while x < floatmax(T)
#         test = parse(T, Ryu.writeshortest(x)) == x
#         if !test

#             fails += 1
#         else
#             success += 1
#         end
#         x = nextfloat(x)
#         i += 1

#     end
#     return fails / (fails + success)
# end

end # Float16

@testset "Ryu.writefixed" begin
    @testset "Basic" begin
        @test Ryu.writefixed(todouble(false, 1234, 99999), 0) ==
            "3291009114715486435425664845573426149758869524108446525879746560"
    end
    @testset "Zero" begin
        @test Ryu.writefixed(0.0, 4) == "0.0000"
        @test Ryu.writefixed(0.0, 3) == "0.000"
        @test Ryu.writefixed(0.0, 2) == "0.00"
        @test Ryu.writefixed(0.0, 1) == "0.0"
        @test Ryu.writefixed(0.0, 0) == "0"
    end
    @testset "MinMax" begin
        @test Ryu.writefixed(todouble(false, 0, 1), 1074) ==
            "0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" *
            "000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" *
            "000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" *
            "000000000000000000000000000000000000000000000000000000049406564584124654417656879286822137" *
            "236505980261432476442558568250067550727020875186529983636163599237979656469544571773092665" *
            "671035593979639877479601078187812630071319031140452784581716784898210368871863605699873072" *
            "305000638740915356498438731247339727316961514003171538539807412623856559117102665855668676" *
            "818703956031062493194527159149245532930545654440112748012970999954193198940908041656332452" *
            "475714786901472678015935523861155013480352649347201937902681071074917033322268447533357208" *
            "324319360923828934583680601060115061698097530783422773183292479049825247307763759272478746" *
            "560847782037344696995336470179726777175851256605511991315048911014510378627381672509558373" *
            "89733598993664809941164205702637090279242767544565229087538682506419718265533447265625"
        @test Ryu.writefixed(todouble(false, 2046, 0xFFFFFFFFFFFFF), 0) ==
            "179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558" *
            "632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245" *
            "490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168" *
            "738177180919299881250404026184124858368"
    end
    @testset "RoundToEven" begin
        @test Ryu.writefixed(0.125, 3) == "0.125"
        @test Ryu.writefixed(0.125, 2) == "0.12"
        @test Ryu.writefixed(0.375, 3) == "0.375"
        @test Ryu.writefixed(0.375, 2) == "0.38"
    end
    @testset "RoundToEvenInteger" begin
        @test Ryu.writefixed(2.5, 1) == "2.5"
        @test Ryu.writefixed(2.5, 0) == "2"
        @test Ryu.writefixed(3.5, 1) == "3.5"
        @test Ryu.writefixed(3.5, 0) == "4"
    end
    @testset "NonRoundToEvenScenarios" begin
        @test Ryu.writefixed(0.748046875, 3) == "0.748"
        @test Ryu.writefixed(0.748046875, 2) == "0.75"
        @test Ryu.writefixed(0.748046875, 1) == "0.7"

        @test Ryu.writefixed(0.2509765625, 3) == "0.251"
        @test Ryu.writefixed(0.2509765625, 2) == "0.25"
        @test Ryu.writefixed(0.2509765625, 1) == "0.3"

        @test Ryu.writefixed(todouble(false, 1021, 1), 54) == "0.250000000000000055511151231257827021181583404541015625"
        @test Ryu.writefixed(todouble(false, 1021, 1),  3) == "0.250"
        @test Ryu.writefixed(todouble(false, 1021, 1),  2) == "0.25"
        @test Ryu.writefixed(todouble(false, 1021, 1),  1) == "0.3"
    end
    @testset "VaryingPrecision" begin
        @test Ryu.writefixed(1729.142857142857, 47) == "1729.14285714285711037518922239542007446289062500000"
        @test Ryu.writefixed(1729.142857142857, 46) == "1729.1428571428571103751892223954200744628906250000"
        @test Ryu.writefixed(1729.142857142857, 45) == "1729.142857142857110375189222395420074462890625000"
        @test Ryu.writefixed(1729.142857142857, 44) == "1729.14285714285711037518922239542007446289062500"
        @test Ryu.writefixed(1729.142857142857, 43) == "1729.1428571428571103751892223954200744628906250"
        @test Ryu.writefixed(1729.142857142857, 42) == "1729.142857142857110375189222395420074462890625"
        @test Ryu.writefixed(1729.142857142857, 41) == "1729.14285714285711037518922239542007446289062"
        @test Ryu.writefixed(1729.142857142857, 40) == "1729.1428571428571103751892223954200744628906"
        @test Ryu.writefixed(1729.142857142857, 39) == "1729.142857142857110375189222395420074462891"
        @test Ryu.writefixed(1729.142857142857, 38) == "1729.14285714285711037518922239542007446289"
        @test Ryu.writefixed(1729.142857142857, 37) == "1729.1428571428571103751892223954200744629"
        @test Ryu.writefixed(1729.142857142857, 36) == "1729.142857142857110375189222395420074463"
        @test Ryu.writefixed(1729.142857142857, 35) == "1729.14285714285711037518922239542007446"
        @test Ryu.writefixed(1729.142857142857, 34) == "1729.1428571428571103751892223954200745"
        @test Ryu.writefixed(1729.142857142857, 33) == "1729.142857142857110375189222395420074"
        @test Ryu.writefixed(1729.142857142857, 32) == "1729.14285714285711037518922239542007"
        @test Ryu.writefixed(1729.142857142857, 31) == "1729.1428571428571103751892223954201"
        @test Ryu.writefixed(1729.142857142857, 30) == "1729.142857142857110375189222395420"
        @test Ryu.writefixed(1729.142857142857, 29) == "1729.14285714285711037518922239542"
        @test Ryu.writefixed(1729.142857142857, 28) == "1729.1428571428571103751892223954"
        @test Ryu.writefixed(1729.142857142857, 27) == "1729.142857142857110375189222395"
        @test Ryu.writefixed(1729.142857142857, 26) == "1729.14285714285711037518922240"
        @test Ryu.writefixed(1729.142857142857, 25) == "1729.1428571428571103751892224"
        @test Ryu.writefixed(1729.142857142857, 24) == "1729.142857142857110375189222"
        @test Ryu.writefixed(1729.142857142857, 23) == "1729.14285714285711037518922"
        @test Ryu.writefixed(1729.142857142857, 22) == "1729.1428571428571103751892"
        @test Ryu.writefixed(1729.142857142857, 21) == "1729.142857142857110375189"
        @test Ryu.writefixed(1729.142857142857, 20) == "1729.14285714285711037519"
        @test Ryu.writefixed(1729.142857142857, 19) == "1729.1428571428571103752"
        @test Ryu.writefixed(1729.142857142857, 18) == "1729.142857142857110375"
        @test Ryu.writefixed(1729.142857142857, 17) == "1729.14285714285711038"
        @test Ryu.writefixed(1729.142857142857, 16) == "1729.1428571428571104"
        @test Ryu.writefixed(1729.142857142857, 15) == "1729.142857142857110"
        @test Ryu.writefixed(1729.142857142857, 14) == "1729.14285714285711"
        @test Ryu.writefixed(1729.142857142857, 13) == "1729.1428571428571"
        @test Ryu.writefixed(1729.142857142857, 12) == "1729.142857142857"
        @test Ryu.writefixed(1729.142857142857, 11) == "1729.14285714286"
        @test Ryu.writefixed(1729.142857142857, 10) == "1729.1428571429"
        @test Ryu.writefixed(1729.142857142857,  9) == "1729.142857143"
        @test Ryu.writefixed(1729.142857142857,  8) == "1729.14285714"
        @test Ryu.writefixed(1729.142857142857,  7) == "1729.1428571"
        @test Ryu.writefixed(1729.142857142857,  6) == "1729.142857"
        @test Ryu.writefixed(1729.142857142857,  5) == "1729.14286"
        @test Ryu.writefixed(1729.142857142857,  4) == "1729.1429"
        @test Ryu.writefixed(1729.142857142857,  3) == "1729.143"
        @test Ryu.writefixed(1729.142857142857,  2) == "1729.14"
        @test Ryu.writefixed(1729.142857142857,  1) == "1729.1"
        @test Ryu.writefixed(1729.142857142857,  0) == "1729"
    end

    @testset "Carrying" begin
        @test Ryu.writefixed(  0.0009, 4) == "0.0009"
        @test Ryu.writefixed(  0.0009, 3) == "0.001"
        @test Ryu.writefixed(  0.0029, 4) == "0.0029"
        @test Ryu.writefixed(  0.0029, 3) == "0.003"
        @test Ryu.writefixed(  0.0099, 4) == "0.0099"
        @test Ryu.writefixed(  0.0099, 3) == "0.010"
        @test Ryu.writefixed(  0.0299, 4) == "0.0299"
        @test Ryu.writefixed(  0.0299, 3) == "0.030"
        @test Ryu.writefixed(  0.0999, 4) == "0.0999"
        @test Ryu.writefixed(  0.0999, 3) == "0.100"
        @test Ryu.writefixed(  0.2999, 4) == "0.2999"
        @test Ryu.writefixed(  0.2999, 3) == "0.300"
        @test Ryu.writefixed(  0.9999, 4) == "0.9999"
        @test Ryu.writefixed(  0.9999, 3) == "1.000"
        @test Ryu.writefixed(  2.9999, 4) == "2.9999"
        @test Ryu.writefixed(  2.9999, 3) == "3.000"
        @test Ryu.writefixed(  9.9999, 4) == "9.9999"
        @test Ryu.writefixed(  9.9999, 3) == "10.000"
        @test Ryu.writefixed( 29.9999, 4) == "29.9999"
        @test Ryu.writefixed( 29.9999, 3) == "30.000"
        @test Ryu.writefixed( 99.9999, 4) == "99.9999"
        @test Ryu.writefixed( 99.9999, 3) == "100.000"
        @test Ryu.writefixed(299.9999, 4) == "299.9999"
        @test Ryu.writefixed(299.9999, 3) == "300.000"

        @test Ryu.writefixed(  0.09, 2) == "0.09"
        @test Ryu.writefixed(  0.09, 1) == "0.1"
        @test Ryu.writefixed(  0.29, 2) == "0.29"
        @test Ryu.writefixed(  0.29, 1) == "0.3"
        @test Ryu.writefixed(  0.99, 2) == "0.99"
        @test Ryu.writefixed(  0.99, 1) == "1.0"
        @test Ryu.writefixed(  2.99, 2) == "2.99"
        @test Ryu.writefixed(  2.99, 1) == "3.0"
        @test Ryu.writefixed(  9.99, 2) == "9.99"
        @test Ryu.writefixed(  9.99, 1) == "10.0"
        @test Ryu.writefixed( 29.99, 2) == "29.99"
        @test Ryu.writefixed( 29.99, 1) == "30.0"
        @test Ryu.writefixed( 99.99, 2) == "99.99"
        @test Ryu.writefixed( 99.99, 1) == "100.0"
        @test Ryu.writefixed(299.99, 2) == "299.99"
        @test Ryu.writefixed(299.99, 1) == "300.0"

        @test Ryu.writefixed(  0.9, 1) == "0.9"
        @test Ryu.writefixed(  0.9, 0) == "1"
        @test Ryu.writefixed(  2.9, 1) == "2.9"
        @test Ryu.writefixed(  2.9, 0) == "3"
        @test Ryu.writefixed(  9.9, 1) == "9.9"
        @test Ryu.writefixed(  9.9, 0) == "10"
        @test Ryu.writefixed( 29.9, 1) == "29.9"
        @test Ryu.writefixed( 29.9, 0) == "30"
        @test Ryu.writefixed( 99.9, 1) == "99.9"
        @test Ryu.writefixed( 99.9, 0) == "100"
        @test Ryu.writefixed(299.9, 1) == "299.9"
        @test Ryu.writefixed(299.9, 0) == "300"
    end

    @testset "RoundingResultZero" begin
        @test Ryu.writefixed(0.004, 3) == "0.004"
        @test Ryu.writefixed(0.004, 2) == "0.00"
        @test Ryu.writefixed(0.4, 1) == "0.4"
        @test Ryu.writefixed(0.4, 0) == "0"
        @test Ryu.writefixed(0.5, 1) == "0.5"
        @test Ryu.writefixed(0.5, 0) == "0"
    end

    @testset "Regression" begin
        @test Ryu.writefixed(7.018232e-82, 6) == "0.000000"
    end

end # fixed

@testset "Ryu.writeexp" begin

@testset "Basic" begin
    @test Ryu.writeexp(todouble(false, 1234, 99999), 62) ==
    "3.29100911471548643542566484557342614975886952410844652587974656e+63"
end

@testset "Zero" begin
    @test Ryu.writeexp(0.0, 4) == "0.0000e+00"
    @test Ryu.writeexp(0.0, 3) == "0.000e+00"
    @test Ryu.writeexp(0.0, 2) == "0.00e+00"
    @test Ryu.writeexp(0.0, 1) == "0.0e+00"
    @test Ryu.writeexp(0.0, 0) == "0e+00"
end

@testset "MinMax" begin
    @test Ryu.writeexp(todouble(false, 0, 1), 750) ==
    "4.9406564584124654417656879286822137236505980261432476442558568250067550727020875186529983" *
    "636163599237979656469544571773092665671035593979639877479601078187812630071319031140452784" *
    "581716784898210368871863605699873072305000638740915356498438731247339727316961514003171538" *
    "539807412623856559117102665855668676818703956031062493194527159149245532930545654440112748" *
    "012970999954193198940908041656332452475714786901472678015935523861155013480352649347201937" *
    "902681071074917033322268447533357208324319360923828934583680601060115061698097530783422773" *
    "183292479049825247307763759272478746560847782037344696995336470179726777175851256605511991" *
    "315048911014510378627381672509558373897335989936648099411642057026370902792427675445652290" *
    "87538682506419718265533447265625e-324"

    @test Ryu.writeexp(todouble(false, 2046, 0xFFFFFFFFFFFFF), 308) ==
    "1.7976931348623157081452742373170435679807056752584499659891747680315726078002853876058955" *
    "863276687817154045895351438246423432132688946418276846754670353751698604991057655128207624" *
    "549009038932894407586850845513394230458323690322294816580855933212334827479782620414472316" *
    "8738177180919299881250404026184124858368e+308"
end

@testset "RoundToEven" begin
    @test Ryu.writeexp(0.125, 2) == "1.25e-01"
    @test Ryu.writeexp(0.125, 1) == "1.2e-01"
    @test Ryu.writeexp(0.375, 2) == "3.75e-01"
    @test Ryu.writeexp(0.375, 1) == "3.8e-01"
end

@testset "RoundToEvenInteger" begin
    @test Ryu.writeexp(2.5, 1) == "2.5e+00"
    @test Ryu.writeexp(2.5, 0) == "2e+00"
    @test Ryu.writeexp(3.5, 1) == "3.5e+00"
    @test Ryu.writeexp(3.5, 0) == "4e+00"
end

@testset "NonRoundToEvenScenarios" begin
    @test Ryu.writeexp(0.748046875, 2) == "7.48e-01"
    @test Ryu.writeexp(0.748046875, 1) == "7.5e-01"
    @test Ryu.writeexp(0.748046875, 0) == "7e-01"    # 0.75 would round to "8e-01", but this is smaller

    @test Ryu.writeexp(0.2509765625, 2) == "2.51e-01"
    @test Ryu.writeexp(0.2509765625, 1) == "2.5e-01"
    @test Ryu.writeexp(0.2509765625, 0) == "3e-01"    # 0.25 would round to "2e-01", but this is larger

    @test Ryu.writeexp(todouble(false, 1021, 1), 53) ==
    "2.50000000000000055511151231257827021181583404541015625e-01"
    @test Ryu.writeexp(todouble(false, 1021, 1),  2) ==
    "2.50e-01"
    @test Ryu.writeexp(todouble(false, 1021, 1),  1) ==
    "2.5e-01"
    @test Ryu.writeexp(todouble(false, 1021, 1),  0) ==
    "3e-01"    # 0.25 would round to "2e-01", but this is larger (again)
end

@testset "VaryingPrecision" begin
    @test Ryu.writeexp(1729.142857142857, 50) == "1.72914285714285711037518922239542007446289062500000e+03"
    @test Ryu.writeexp(1729.142857142857, 49) == "1.7291428571428571103751892223954200744628906250000e+03"
    @test Ryu.writeexp(1729.142857142857, 48) == "1.729142857142857110375189222395420074462890625000e+03"
    @test Ryu.writeexp(1729.142857142857, 47) == "1.72914285714285711037518922239542007446289062500e+03"
    @test Ryu.writeexp(1729.142857142857, 46) == "1.7291428571428571103751892223954200744628906250e+03"
    @test Ryu.writeexp(1729.142857142857, 45) == "1.729142857142857110375189222395420074462890625e+03"
    @test Ryu.writeexp(1729.142857142857, 44) == "1.72914285714285711037518922239542007446289062e+03"
    @test Ryu.writeexp(1729.142857142857, 43) == "1.7291428571428571103751892223954200744628906e+03"
    @test Ryu.writeexp(1729.142857142857, 42) == "1.729142857142857110375189222395420074462891e+03"
    @test Ryu.writeexp(1729.142857142857, 41) == "1.72914285714285711037518922239542007446289e+03"
    @test Ryu.writeexp(1729.142857142857, 40) == "1.7291428571428571103751892223954200744629e+03"
    @test Ryu.writeexp(1729.142857142857, 39) == "1.729142857142857110375189222395420074463e+03"
    @test Ryu.writeexp(1729.142857142857, 38) == "1.72914285714285711037518922239542007446e+03"
    @test Ryu.writeexp(1729.142857142857, 37) == "1.7291428571428571103751892223954200745e+03"
    @test Ryu.writeexp(1729.142857142857, 36) == "1.729142857142857110375189222395420074e+03"
    @test Ryu.writeexp(1729.142857142857, 35) == "1.72914285714285711037518922239542007e+03"
    @test Ryu.writeexp(1729.142857142857, 34) == "1.7291428571428571103751892223954201e+03"
    @test Ryu.writeexp(1729.142857142857, 33) == "1.729142857142857110375189222395420e+03"
    @test Ryu.writeexp(1729.142857142857, 32) == "1.72914285714285711037518922239542e+03"
    @test Ryu.writeexp(1729.142857142857, 31) == "1.7291428571428571103751892223954e+03"
    @test Ryu.writeexp(1729.142857142857, 30) == "1.729142857142857110375189222395e+03"
    @test Ryu.writeexp(1729.142857142857, 29) == "1.72914285714285711037518922240e+03"
    @test Ryu.writeexp(1729.142857142857, 28) == "1.7291428571428571103751892224e+03"
    @test Ryu.writeexp(1729.142857142857, 27) == "1.729142857142857110375189222e+03"
    @test Ryu.writeexp(1729.142857142857, 26) == "1.72914285714285711037518922e+03"
    @test Ryu.writeexp(1729.142857142857, 25) == "1.7291428571428571103751892e+03"
    @test Ryu.writeexp(1729.142857142857, 24) == "1.729142857142857110375189e+03"
    @test Ryu.writeexp(1729.142857142857, 23) == "1.72914285714285711037519e+03"
    @test Ryu.writeexp(1729.142857142857, 22) == "1.7291428571428571103752e+03"
    @test Ryu.writeexp(1729.142857142857, 21) == "1.729142857142857110375e+03"
    @test Ryu.writeexp(1729.142857142857, 20) == "1.72914285714285711038e+03"
    @test Ryu.writeexp(1729.142857142857, 19) == "1.7291428571428571104e+03"
    @test Ryu.writeexp(1729.142857142857, 18) == "1.729142857142857110e+03"
    @test Ryu.writeexp(1729.142857142857, 17) == "1.72914285714285711e+03"
    @test Ryu.writeexp(1729.142857142857, 16) == "1.7291428571428571e+03"
    @test Ryu.writeexp(1729.142857142857, 15) == "1.729142857142857e+03"
    @test Ryu.writeexp(1729.142857142857, 14) == "1.72914285714286e+03"
    @test Ryu.writeexp(1729.142857142857, 13) == "1.7291428571429e+03"
    @test Ryu.writeexp(1729.142857142857, 12) == "1.729142857143e+03"
    @test Ryu.writeexp(1729.142857142857, 11) == "1.72914285714e+03"
    @test Ryu.writeexp(1729.142857142857, 10) == "1.7291428571e+03"
    @test Ryu.writeexp(1729.142857142857,  9) == "1.729142857e+03"
    @test Ryu.writeexp(1729.142857142857,  8) == "1.72914286e+03"
    @test Ryu.writeexp(1729.142857142857,  7) == "1.7291429e+03"
    @test Ryu.writeexp(1729.142857142857,  6) == "1.729143e+03"
    @test Ryu.writeexp(1729.142857142857,  5) == "1.72914e+03"
    @test Ryu.writeexp(1729.142857142857,  4) == "1.7291e+03"
    @test Ryu.writeexp(1729.142857142857,  3) == "1.729e+03"
    @test Ryu.writeexp(1729.142857142857,  2) == "1.73e+03"
    @test Ryu.writeexp(1729.142857142857,  1) == "1.7e+03"
    @test Ryu.writeexp(1729.142857142857,  0) == "2e+03"
end

@testset "Carrying" begin
    @test Ryu.writeexp(2.0009, 4) == "2.0009e+00"
    @test Ryu.writeexp(2.0009, 3) == "2.001e+00"
    @test Ryu.writeexp(2.0029, 4) == "2.0029e+00"
    @test Ryu.writeexp(2.0029, 3) == "2.003e+00"
    @test Ryu.writeexp(2.0099, 4) == "2.0099e+00"
    @test Ryu.writeexp(2.0099, 3) == "2.010e+00"
    @test Ryu.writeexp(2.0299, 4) == "2.0299e+00"
    @test Ryu.writeexp(2.0299, 3) == "2.030e+00"
    @test Ryu.writeexp(2.0999, 4) == "2.0999e+00"
    @test Ryu.writeexp(2.0999, 3) == "2.100e+00"
    @test Ryu.writeexp(2.2999, 4) == "2.2999e+00"
    @test Ryu.writeexp(2.2999, 3) == "2.300e+00"
    @test Ryu.writeexp(2.9999, 4) == "2.9999e+00"
    @test Ryu.writeexp(2.9999, 3) == "3.000e+00"
    @test Ryu.writeexp(9.9999, 4) == "9.9999e+00"
    @test Ryu.writeexp(9.9999, 3) == "1.000e+01"

    @test Ryu.writeexp(2.09, 2) == "2.09e+00"
    @test Ryu.writeexp(2.09, 1) == "2.1e+00"
    @test Ryu.writeexp(2.29, 2) == "2.29e+00"
    @test Ryu.writeexp(2.29, 1) == "2.3e+00"
    @test Ryu.writeexp(2.99, 2) == "2.99e+00"
    @test Ryu.writeexp(2.99, 1) == "3.0e+00"
    @test Ryu.writeexp(9.99, 2) == "9.99e+00"
    @test Ryu.writeexp(9.99, 1) == "1.0e+01"

    @test Ryu.writeexp(2.9, 1) == "2.9e+00"
    @test Ryu.writeexp(2.9, 0) == "3e+00"
    @test Ryu.writeexp(9.9, 1) == "9.9e+00"
    @test Ryu.writeexp(9.9, 0) == "1e+01"
end

@testset "Exponents" begin
    @test Ryu.writeexp(9.99e-100, 2) == "9.99e-100"
    @test Ryu.writeexp(9.99e-99 , 2) == "9.99e-99"
    @test Ryu.writeexp(9.99e-10 , 2) == "9.99e-10"
    @test Ryu.writeexp(9.99e-09 , 2) == "9.99e-09"
    @test Ryu.writeexp(9.99e-01 , 2) == "9.99e-01"
    @test Ryu.writeexp(9.99e+00 , 2) == "9.99e+00"
    @test Ryu.writeexp(9.99e+01 , 2) == "9.99e+01"
    @test Ryu.writeexp(9.99e+09 , 2) == "9.99e+09"
    @test Ryu.writeexp(9.99e+10 , 2) == "9.99e+10"
    @test Ryu.writeexp(9.99e+99 , 2) == "9.99e+99"
    @test Ryu.writeexp(9.99e+100, 2) == "9.99e+100"

    @test Ryu.writeexp(9.99e-100, 1) == "1.0e-99"
    @test Ryu.writeexp(9.99e-99 , 1) == "1.0e-98"
    @test Ryu.writeexp(9.99e-10 , 1) == "1.0e-09"
    @test Ryu.writeexp(9.99e-09 , 1) == "1.0e-08"
    @test Ryu.writeexp(9.99e-01 , 1) == "1.0e+00"
    @test Ryu.writeexp(9.99e+00 , 1) == "1.0e+01"
    @test Ryu.writeexp(9.99e+01 , 1) == "1.0e+02"
    @test Ryu.writeexp(9.99e+09 , 1) == "1.0e+10"
    @test Ryu.writeexp(9.99e+10 , 1) == "1.0e+11"
    @test Ryu.writeexp(9.99e+99 , 1) == "1.0e+100"
    @test Ryu.writeexp(9.99e+100, 1) == "1.0e+101"
end

@testset "PrintDecimalPoint" begin
  # These values exercise each codepath.
    @test Ryu.writeexp(1e+54, 0) == "1e+54"
    @test Ryu.writeexp(1e+54, 1) == "1.0e+54"
    @test Ryu.writeexp(1e-63, 0) == "1e-63"
    @test Ryu.writeexp(1e-63, 1) == "1.0e-63"
    @test Ryu.writeexp(1e+83, 0) == "1e+83"
    @test Ryu.writeexp(1e+83, 1) == "1.0e+83"
end

end # exp

end # Ryu

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: milesfrain/julia

max_forks_repo_name: milesfrain/julia


text:
# This file is a part of MJDSigGen, licensed under the MIT License (MIT).

__precompile__(true)

module MJDSigGen

# Make Pkg aware that deps.jl will require Libdl:
import Libdl

const deps_jl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps_jl)
    include(deps_jl)
else
    error("MJDSigGen is not properly installed. Run Pkg.build(\"MJDSigGen\").")
end

macro sgsym(func)
    :(($(esc(func)), libmjd_siggen))
end

import Base

export SigGenSetup, read_config!, read_config, fieldgen, read_fields, signal_calc_init!,
    signal_calc_init, get_signal!, outside_detector, with_coll_effects!, coll_effects_off!,
    getδτ, getδτ!

using DelimitedFiles: readdlm

include("types.jl")
include("util.jl")
include("wrappers.jl")

end # module

ext: jl


text:
#=
   BLOG example/uniform-ball.blog
   """
   Model file for balls in an urn, allowing observation errors.
   This version uses a Poisson prior for the number of balls.
   """

   Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

      numBalls    2.8192    1.0196     0.0102    0.0225   1557.7583    0.9999      464.5864
      ...
             q    0.3318    0.4709     0.0047    0.0114   1815.4962    1.0000      541.4543
      obsPost1    0.6332    0.4820     0.0048    0.0074   5313.9656    1.0000     1584.8391

   Distributions of variable numBalls (num:0)
   4.00000 =>    3325  (0.332500)
   2.00000 =>    2869  (0.286900)
   3.00000 =>    2674  (0.267400)
   1.00000 =>    1132  (0.113200)

  Cf ~/webppl/uniform_ball.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function uniform_ball(obs=[1,0,1,0,1,0,1,0])
    n = length(obs)

    numBalls ~ DiscreteUniform(1,4)

    isBlue ~ filldist(flip(0.5),n)

    ballDrawn ~ filldist(DiscreteUniform(1,numBalls),n)

    for d in 1:n
        if isBlue[ballDrawn[d]]
            obs[d] ~ flip(0.8)
        else
            obs[d] ~ flip(0.2)
        end
    end

    q ~ Dirac(ballDrawn[1] == ballDrawn[2])

    if isBlue[ballDrawn[1]]
        obsPost1 ~ flip(0.8)
    else
        obsPost1 ~ flip(0.2)
    end

end

model = uniform_ball()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:numBalls)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Reagular print
print("Hello world\n")

# a println statement
println("Hello world - line!")

# Prining a variable value
x = 42
println("The value of x = ", x)

# printing multiple values
y = 11
println("The value of x = ", x, ", and the value of y = ", y)

# printing a long line
print("The value of x = ")
print(x)
print(", And the value of y = ")
println(y)

# formatting a float value
println("pi = ", Float64(pi))
using Printf
@printf("%s %.3f\n", "pi to 3 decimal places = ", Float64(pi))

#concatenating string variables
x_str = "My name is"
y_str = "Zizzybaluba"
sentence = x_str * " " * y_str
println(sentence)

# we ca also concatentate using string
sentence2 = string(x_str, " ", y_str)
println("this is sentence2: ", sentence2)


text:
using RCall

#Run the KD simulations in parallel
#change the number of cores used as needed
using Distributed
addprocs(4)
#The @everywhere macro runs the command/defines the thing on all the processors
@everywhere include("ProblemGenerator.jl")
@everywhere include("VirusCallBacks.jl")

########################################################################
#Knockdown simulations of IRF7 and TREX1 (Figure 4)
########################################################################

#Putting the calculations in a function will allow for easy parallelization
@everywhere function KnockDownSim(par2Change,prob)
  #Percentages to run the KD simulations
  KnockDownVals=1:-0.25:0
  percentLabels = [string(i)*"%" for i=0:25:100]
  kdSamples = length(KnockDownVals)

  #Vector to save all IFN dynamics from simulation
  solKD =  Vector(undef,kdSamples)

  #Loop through the KD percents and solve ODEs
  for (i,percent) in enumerate(KnockDownVals)
      #Copy the parameters and change desired one
      θcurrent = deepcopy(prob.p)
      #Check if running ABM or ODE
      if θcurrent isa Array #Are the parameters from the ODE model?
          θcurrent[par2Change] = θcurrent[par2Change] * percent
          #Print to track progress
          println(θcurrent[[20,23]])
      else #Oh nope its the ABM model then
          θcurrent.par[par2Change] = θcurrent.par[par2Change] * percent
          #Print to track progress
          println(θcurrent.par[[20,23]])
      end

      #Redefine problem with the new parameters and solve
      probKD = remake(prob; p=θcurrent)

      if θcurrent isa Array
        sol = solve(probKD,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)
        solKD[i] = vec(sol[7,:])
      else
        #Check if ISD or Virus infection
        if θcurrent.DNAReplicate == 0 #ISD
          sol = solve(probKD,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)
        else #Virus
          sol = solve(probKD,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback = cb)
        end
        solKD[i] =vec(sol[:,:,7,:])
      end

    end

  #Return the IFN dynamics
  return solKD
end

########################################################################
#Knockdown ABM simulations (Deterministic and Homogeneous)
########################################################################

#Define a problem
@everywhere prob = ModelSetup(:ISD,:notStochastic,:Homo)
#Run the KD simulation
KDSols = pmap(x-> KnockDownSim(x,prob),[20,23])

#Save the simulation in a Dataframe to pass on to R (for plotting)
KDDataDet = DataFrame()

#Seperate out the columns for each KD
KDDataDet.IRF7KD = vcat(KDSols[1]...)
KDDataDet.TREXKD = vcat(KDSols[2]...)

#Create columns that keep track of cell ID, time, and KD percent
timeLength = 481
KnockDownVals=1:-0.25:0
percentLabels = [string(i)*"%" for i=0:25:100]
kdSamples = length(KnockDownVals)

#Put everything into the DataFrame
KDDataDet.Cell = repeat(1:nCells,timeLength*kdSamples)
KDDataDet.Time = repeat(0:0.1:48,inner=nCells,outer=kdSamples)
KDDataDet.Percent = repeat(KnockDownVals,inner=nCells*timeLength)

########################################################################
#Knockdown ABM simulations (Stochastic and Heterogeneous)
########################################################################
#Define a problem
@everywhere prob = ModelSetup(:ISD,:Stochastic,:Hetero)
#Run the KD simulation
KDSols = pmap(x-> KnockDownSim(x,prob),[20,23])

#Save the simulation in a Dataframe to pass on to R (for plotting)
KDDataStoch = DataFrame()

#Seperate out the columns for each KD
KDDataStoch.IRF7KD = vcat(KDSols[1]...)
KDDataStoch.TREXKD = vcat(KDSols[2]...)

#Put everything into the DataFrame
KDDataStoch.Cell = repeat(1:nCells,timeLength*kdSamples)
KDDataStoch.Time = repeat(0:0.1:48,inner=nCells,outer=kdSamples)
KDDataStoch.Percent = repeat(KnockDownVals,inner=nCells*timeLength)

########################################################################
#Knockdown ODE simulations
########################################################################

#ODE IC
ODEu0 = zeros(species)
ODEu0[[1,2,3,5]] .= m2c([1e3,1e3,1e3, 1e4])

#ODE Pars
ODEpars = [2.6899, 4.8505, 0.0356, 7.487, 517.4056, 22328.3852, 11226.3682,0.9341,
         206.9446, 10305.461, 47639.70295,3.8474, 13.006, 78.2048, 0.0209,
         0.0059, 0.001, 0.0112, 0.001, 99.9466, 15.1436,0.0276, 237539.3249,
         61688.259, 0.96, 0.347, 12.2428736,1.2399, 1.5101, 0.347, 0.165, 6.9295,
         0.0178]
 #Append  cGAStot, Stingtot, IRF3tot to the parameters
 append!(ODEpars,ODEu0[[1,3,5]])

 #ODE Model
 function ODEmodel(du,u,p,t)
   #Species
   cGAS, DNA, Sting, cGAMP, IRF3, IFNβm, IFNβ, STAT, SOCSm, IRF7m, TREX1m, IRF7, TREX1 = u

   #Parameters
   k1f, k1r, k3f, k3r, k4f, kcat5, Km5, k5r, kcat6, Km6, kcat7, Km7, kcat8, Km8, k8f, k9f, k10f1, k10f2, k11f, k12f, k13f, k6f, kcat2, Km2, τ4, τ6, τ7, τ8, τ9, τ10, τ11, τ12, τ13, cGAStot, Stingtot, IRF3tot = p

   #Update derivatives for each species according to model
   du[1] = -k1f*cGAS*DNA + k1r*(cGAStot - cGAS)
   du[2] = -k1f*cGAS*DNA + k1r*(cGAStot - cGAS) - kcat2*TREX1*DNA / (Km2 + DNA)
   du[3] = -k3f*cGAMP*Sting + k3r*(Stingtot - Sting)
   du[4] = k4f*(cGAStot - cGAS) - k3f*cGAMP*Sting + k3f*(Stingtot - Sting) - τ4*cGAMP
   du[5] = -kcat5*IRF3*(Stingtot - Sting) / (Km5 +IRF3) + k5r*(IRF3tot - IRF3)
   du[6] = kcat6*(IRF3tot - IRF3) / (Km6 + (IRF3tot - IRF3)) + k6f*IRF7 - τ6*IFNβm
   du[7] = kcat7*IFNβm / (Km7 + IFNβm) - τ7*IFNβ
   du[8] = kcat8*IFNβ / (Km8 + IFNβ) * 1.0/(1.0+k8f*SOCSm) - τ8*STAT
   du[9] = k9f*STAT - τ9*SOCSm
   du[10] = k10f1*STAT + k10f2*IRF7 - τ10*IRF7m
   du[11] = k11f*STAT - τ11*TREX1m
   du[12] = k12f*IRF7m - τ12*IRF7
   du[13] = k13f*TREX1m - τ13*TREX1
 end

 ODEprob = ODEProblem(ODEmodel,ODEu0,tspan,ODEpars)

 #Run the KD simulation (THis could also be parallelized, but its really fast)
 ODEsols = map(x-> KnockDownSim(x,ODEprob),[20,23])

#GEnerate the DataFrame for the ODE results
 KDDataODE = DataFrame()
 KDDataODE.IRF7KD = vcat(ODEsols[1]...)
 KDDataODE.TREXKD = vcat(ODEsols[2]...)
 KDDataODE.Time = repeat(0:0.1:48,outer=kdSamples)
 KDDataODE.Percent = repeat(KnockDownVals,inner=timeLength)

 ########################################################################
 #Plotting the ABM and ODE solutions
 ########################################################################

 @rput KDDataDet
 @rput KDDataStoch
 @rput KDDataODE
 R"""
 library(ggplot2)
 library(ggpubr)

#--------------Det--------------
 ISDDetAve = aggregate(KDDataDet[,1:2], list(KDDataDet$Percent,KDDataDet$Time), mean)
 colnames(ISDDetAve)[1] <- "Percent"
 colnames(ISDDetAve)[2] <- "Time"

 stateSD = aggregate(KDDataDet[,1:2], list(KDDataDet$Percent,KDDataDet$Time), sd)
  lowDet = ISDDetAve[,3:4] -  stateSD[,3:4]
  highDet = ISDDetAve[,3:4] + stateSD[,3:4]
 #low = aggregate(KDDataDet[,1:2], list(KDDataDet$Percent,KDDataDet$Time), FUN = 'quantile',probs=0.05)
 #high = aggregate(KDDataDet[,1:2], list(KDDataDet$Percent,KDDataDet$Time), FUN = 'quantile',probs=0.95)

#--------------Stoch--------------
 ISDStochAve = aggregate(KDDataStoch[,1:2], list(KDDataStoch$Percent,KDDataStoch$Time), mean)
 colnames(ISDStochAve)[1] <- "Percent"
 colnames(ISDStochAve)[2] <- "Time"

 stateSD = aggregate(KDDataStoch[,1:2], list(KDDataStoch$Percent,KDDataStoch$Time), sd)
  lowStoch = ISDStochAve[,3:4] -  stateSD[,3:4]
  highStoch = ISDStochAve[,3:4] + stateSD[,3:4]

 commonFigureOptions <- list(scale_x_continuous(breaks=seq(0, 48, 12)),
   theme_pubr(border=TRUE),
   xlab("Time (hours)"),
   theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1))

 p1 <- ggplot(KDDataODE) + geom_line(aes(y=TREXKD, x=Time, group=factor(Percent), color = factor(Percent))) +
   ggtitle("ODE Model: \n TREX1 Knockdown") +
   ylab("IFN (nM)") +
   commonFigureOptions

 p2 <- ggplot(KDDataODE) + geom_line(aes(y=IRF7KD, x=Time, group=factor(Percent), color = factor(Percent))) +
   ggtitle("ODE Model: \n IRF7 Knockdown") +
   ylab("IFN (nM)") +
   commonFigureOptions

 p3 <- ggplot(ISDDetAve) + geom_line(aes(y=TREXKD, x=Time, group=factor(Percent), color = factor(Percent))) +
   geom_ribbon(aes(ymin=lowDet$TREXKD, ymax=highDet$TREXKD, x=Time,group=factor(Percent),fill = factor(Percent)), alpha = 0.2) +
   ggtitle("ABM: Deterministic + Homogeneous \n TREX1 Knockdown") +
   ylab("Average IFN (nM)") +
   commonFigureOptions

 p4 <- ggplot(ISDDetAve) + geom_line(aes(y=IRF7KD, x=Time, group=factor(Percent), color = factor(Percent))) +
   geom_ribbon(aes(ymin=lowDet$IRF7KD, ymax=highDet$IRF7KD, x=Time,group=factor(Percent),fill = factor(Percent)), alpha = 0.2) +
   ggtitle("ABM: Deterministic + Homogeneous \n IRF7 Knockdown") +
   ylab("Average IFN (nM)") +
   commonFigureOptions

 p5 <- ggplot(ISDStochAve) + geom_line(aes(y=TREXKD, x=Time, group=factor(Percent), color = factor(Percent))) +
   geom_ribbon(aes(ymin=lowStoch$TREXKD, ymax=highStoch$TREXKD, x=Time,group=factor(Percent),fill = factor(Percent)), alpha = 0.2) +
   ggtitle("ABM: Stochastic + Heterogeneous \n TREX1 Knockdown") +
   ylab("Average IFN (nM)") +
   commonFigureOptions

 p6 <- ggplot(ISDStochAve) + geom_line(aes(y=IRF7KD, x=Time, group=factor(Percent), color = factor(Percent))) +
   geom_ribbon(aes(ymin=lowStoch$IRF7KD, ymax=highStoch$IRF7KD, x=Time,group=factor(Percent),fill = factor(Percent)), alpha = 0.2) +
   ggtitle("ABM: Stochastic + Heterogeneous \n IRF7 Knockdown") +
   ylab("Average IFN (nM)") +
   commonFigureOptions

   figure <- ggarrange(p1, p3, p5, p2, p4, p6,
                       labels = "AUTO",
                       common.legend = TRUE, legend = "right",
                       align = "hv",
                       ncol = 3, nrow = 2)

 ggsave("./Figures/Figure4.pdf",width = 15,height=8)
 """


text:
abstract type AbstractPosteriors end

#holds the posterior distributions
mutable struct Posteriors <: AbstractPosteriors
  σ²G!::Function
  ζ²G!::Function
  ζ²P!::Function
  SG!::Function
  SGP!::Function
end

#NOTE: Each posterior distribution as a corresponding test version

function Posteriors(;
  σ²G!::Function = posteriorσ²G!,
  ζ²G!::Function = posteriorζ²G!,
  ζ²P!::Function = posteriorζ²P!,
  SG!::Function = posteriorSG!,
  SGP!::Function = posteriorSGP!)::Posteriors

  return Posteriors(σ²G!, ζ²G!, ζ²P!, SG!, SGP!)
end

function PosteriorsTest(;
  σ²G::Function = posteriorσ²Gtest!,
  ζ²G::Function = posteriorζ²Gtest!,
  ζ²P::Function = posteriorζ²Ptest!,
  SG::Function = posteriorSGtest!,
  SGP::Function = posteriorSGPtest!)::Posteriors

  return Posteriors(σ²G, ζ²G, ζ²P, SG, SGP)
end

#############
function posteriorσ²G!(Θ::Parameters, Π::Priors)::Nothing

  ζ²Gstar::Float64 = 1.0 / (1.0/Θ.ζ²G + 1.0/Θ.ζ²P + 1.0/Π.δ²G)

  μ::Float64 = (Θ.SG/Θ.ζ²G + Θ.SGP/Θ.ζ²P + Π.θG/Π.δ²G) * ζ²Gstar

  Θ.σ²G = rand(TruncatedNormal(μ, ζ²Gstar, 0.0, ∞))

  return nothing
end

#################
function posteriorSG!(Θ::Parameters, Π::Priors)::Nothing
  ζ²SGstar::Float64 = 1.0 / (1.0/Θ.ζ²G + 1.0/Π.δ²SG)

  μ::Float64 = (Θ.σ²G/Θ.ζ²G + Π.θSG/Π.δ²SG) * ζ²SGstar

  Θ.SG = rand(TruncatedNormal(μ, ζ²SGstar, 0.0, ∞))

  return nothing
end

###########################
function posteriorSGP!(Θ::Parameters, Π::Priors)::Nothing
  ζ²SGPstar::Float64 = 1.0 / (1.0/Θ.ζ²G + 1.0/Π.δ²SGP)

  μ::Float64 = (Θ.σ²G/Θ.ζ²G + Π.θSGP/Π.δ²SGP) * ζ²SGPstar

  #Θ.SGP = rand(TruncatedNormal(μ, ζ²SGPstar, 0.0, ∞))
  Θ.SGP = rand(Normal(μ, ζ²SGPstar))

  return nothing
end

#########################
function posteriorζ²G!(Θ::Parameters, Π::Priors)::Nothing
  α::Float64 = Π.αG + 0.5
  β::Float64 = Π.βG + (Θ.SG-Θ.σ²G)^2/2

  Θ.ζ²G = rand(InverseGamma(α, β))

  return nothing
end

#########################
function posteriorζ²P!(Θ::Parameters, Π::Priors)::Nothing
  α::Float64 = Π.αP + 0.5
  β::Float64 = Π.βP + (Θ.SGP-Θ.σ²G)^2/2

  Θ.ζ²P = rand(InverseGamma(α, β))

  return nothing
end

##############****************************#################################
##############****************************#################################
##############****************************#################################
##############****************************#################################

function posteriorσ²Gtest!(Θ::Parameters, Π::Priors)::Nothing

  Θ.σ²G = rand(Normal((Θ.SG/Θ.ζ²G+Θ.SGP/Θ.ζ²P+Π.θG/Π.δ²G)/(1/Θ.ζ²G+1/Θ.ζ²P+1/Π.δ²G),
    1.0/(1/Θ.ζ²G+1/Θ.ζ²P+1/Π.δ²G)))

  return nothing
end

function posteriorSGtest!(Θ::Parameters, Π::Priors)::Nothing

  Θ.SG = rand(Normal((Θ.σ²G/Θ.ζ²G+Π.θSG/Π.δ²SG)/(1/Θ.ζ²G+1/Π.δ²SG),
    1.0/(1/Θ.ζ²G+1/Π.δ²SG)))

  return nothing
end

function posteriorSGPtest!(Θ::Parameters, Π::Priors)::Nothing

  Θ.SGP = rand(Normal((Θ.σ²G/Θ.ζ²G+Π.θSGP/Π.δ²SGP)/(1/Θ.ζ²G+1/Π.δ²SGP),
    1.0/(1/Θ.ζ²G+1/Π.δ²SGP)))

  return nothing
end

function posteriorζ²Gtest!(Θ::Parameters, Π::Priors)::Nothing

  Θ.ζ²G = rand(InverseGamma(Π.αG+0.5, Π.βG + (Θ.SG - Θ.σ²G)^2/2))

  return nothing
end

function posteriorζ²Ptest!(Θ::Parameters, Π::Priors)::Nothing

  Θ.ζ²P = rand(InverseGamma(Π.αP + 0.5, Π.βP + (Θ.SGP - Θ.σ²G)^2/2))

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function partitionnementBiObjectif(instance)
    V, A, P = AEPbiFiltrage(instance)
    graphiqueBis = plt.figure()
    plotPointsNonDomines(V, A, graphiqueBis)
    choix = choixBinaire("\n --> Souhaitez-vous calculer les points non-dominés supportés (o/n) ? ")
    if choix == "o"
        calculPointsSupportes(V, A, graphiqueBis)
    end
    solution = choixSolution(V, A)
    close(graphiqueBis)
    backtracking(instance, V, A, P, solution)
end

function backtracking(instance, V, A, P, solution)
    operationMinIJ = Array{Int64,2}(undef, length(instance.solution.ordreVisiteCamion), length(instance.solution.ordreVisiteCamion))
    for i in 1:length(instance.solution.ordreVisiteCamion)
        for j in 1:length(instance.solution.ordreVisiteCamion)
            operationMinIJ[i,j] = -100
        end
    end

    temp = []
    pred = P[end][solution]
    push!(temp, length(instance.solution.ordreVisiteCamion))
    push!(temp, pred[1])
    operationMinIJ[temp[end], temp[end-1]] = pred[3]
    while pred[1] != 1
        pred = P[pred[1]][pred[2]]
        push!(temp, pred[1])
        operationMinIJ[temp[end], temp[end-1]] = pred[3]
    end
    temp = reverse(temp)

    ordreVisiteCamion = []
    operationsDrone = []
    push!(ordreVisiteCamion, 1)
    for i in 1:length(temp)-1
        operation = operationMinIJ[temp[i],temp[i+1]]
        if operation == -1 # Trajet direct
            noeudCombine = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(ordreVisiteCamion, noeudCombine)
        else # Opération drone
            noeudDrone = instance.solution.ordreVisiteCamion[operation]
            noeudDepart = instance.solution.ordreVisiteCamion[temp[i]]
            noeudFin = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(operationsDrone, [noeudDepart, noeudDrone, noeudFin])
            for l in temp[i]+1:temp[i+1]-1
                noeudActu = instance.solution.ordreVisiteCamion[l]
                if noeudActu != noeudDrone # Noeud camion
                    push!(ordreVisiteCamion, noeudActu)
                end
            end
            push!(ordreVisiteCamion, noeudFin)
        end
    end

    instance.solution.tempsParcours = V[end][solution]
    instance.solution.tempsAttente = A[end][solution]
    instance.solution.ordreVisiteCamion = ordreVisiteCamion
    instance.solution.operationsDrone = operationsDrone
    instance.codeResolution = 5
end

function choixSolution(V, A) # On choisit la solution à afficher (avec laquelle on effectue le backtracking)
    println("\n Liste des points non-dominés :")
    println(" ------------------------------")
    for i in 1:length(V[end])
        println(" --> $i) ($(V[end][i]),$(A[end][i]))")
        text(V[end][i]+1, A[end][i]+1, i)
    end
    choix = choixListe(length(V[end]))
    return parse(Int64, choix)
end

function plotPointsNonDomines(V, A, graphiqueBis)
    ax = graphiqueBis.gca()
    ax.plot(V[end], A[end], marker="x", linestyle="", color="C0", label="Non-dominated points")
    title(L"#$Y_N =$ "*string(length(V[end])))
    grid(true)
    xlabel(L"$z_1$")
    ylabel(L"$z_2$")
    legend(loc="best")
end

function calculPointsSupportes(V, A, graphiqueBis)
    ax = graphiqueBis.gca()
    listeDesPoints = [tuple(V[end][i],A[end][i]) for i in 1:length(V[end])]
    hull = convex_hull_graham_scan!(listeDesPoints) # On doit enlever les points non-supportés qui sont sur les bords de l'enveloppe convexe

    # On cherche les deux points extrêmes
    maxY = -Inf32
    maxIndice = -1
    minY = Inf32
    minIndice = -1
    for i in 1:length(hull)
        if hull[i][2] < minY
            minY = hull[i][2]
            minIndice = i
        end
        if hull[i][2] > maxY
            maxY = hull[i][2]
            maxIndice = i
        end
    end

    # On cherche les points qui sont au-dessus (et donc à enlever)
    vecteurMinMax = [hull[minIndice][1] - hull[maxIndice][1], hull[minIndice][2] - hull[maxIndice][2]]
    hullNew = [hull[minIndice]]
    for i in 1:length(hull)
        if i != minIndice && i != maxIndice
            vecteur = [hull[i][1] - hull[maxIndice][1], hull[i][2] - hull[maxIndice][2]]
            produitVectoriel = (vecteurMinMax[1]*vecteur[2])-(vecteurMinMax[2]*vecteur[1])
            if produitVectoriel < 0 # ok : point supporté
                push!(hullNew, hull[i])
            end
        end
    end
    push!(hullNew, hull[maxIndice])
    hull = hullNew

    YSN = length(hull)
    ax.plot(hull[1][1], hull[1][2], marker="x", color="red", linestyle="", label="Supported non-dominated points")
    for i in 1:length(hull)
        #println(cvh_graham[i])
        ax.plot(hull[i][1], hull[i][2], marker="x", color="red", linestyle="")
        if i != length(hull)
            ax.plot([hull[i][1],hull[i+1][1]], [hull[i][2],hull[i+1][2]], color="red", linestyle="--", linewidth=0.7)
        end
    end

    YN = length(V[end])
    YSN = length(hull)
    YNN = YN - YSN
    title(L"#$Y_{N} =$ "*string(YN)*L"; #$Y_{SN} =$ "*string(YSN)*L"; #$Y_{NN} =$ "*string(YNN))
    legend(loc="best")
end

#################### Bi-objectif ####################

function tempsOperationBi(instance, i,j,k)
    solutionTSP = instance.solution.ordreVisiteCamion # Solution TSP : ordre de visite  du camion
    # Debug
    #println(" Calcul du temps de l'opération : (i=$i, j=$j, k=$k)")
    #println(" Solution TSP : R_tsp = $solutionTSP")
    #print(" r_$i = $(solutionTSP[i]);")
    #print(" r_$j = $(solutionTSP[j]);")
    #if k != -1
    #    println(" r_$k = $(solutionTSP[k])")
    #end
    # Calcul distance parcourue par le drone
    if k == -1 # Pas de noeud drone
        distanceDrone = 0
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        distanceCamion = instance.D[ri,rj]
        tempsAttente = 0.0
    else
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        rk = solutionTSP[k]
        distanceDrone = instance.D[ri,rk] + instance.D[rk,rj]
        distanceCamion = 0
        for l in i:(k-2)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        rkAvant = solutionTSP[k-1]
        rkApres = solutionTSP[k+1]
        distanceCamion += instance.D[rkAvant,rkApres]
        for l in (k+1):(j-1)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        tempsAttente = abs((distanceCamion/instance.vitesseCamion) - (distanceDrone/instance.vitesseDrone))
    end
    tempsParcours = max(distanceCamion/instance.vitesseCamion, distanceDrone/instance.vitesseDrone)
    return tempsParcours, tempsAttente
    #println(" Distance parcourue par le drone : $distanceDrone")
    #println(" Distance parcourue par le camion : $distanceCamion")
end

function calculMeilleuresOperationsBi(instance)
    solutionTSP = instance.solution.ordreVisiteCamion
    tempsMinIJ = Array{Vector{Float64},2}(undef, length(solutionTSP), length(solutionTSP))
    attenteIJ = Array{Vector{Float64},2}(undef, length(solutionTSP), length(solutionTSP))
    operationMinIJ = Array{Vector{Int64},2}(undef, length(solutionTSP), length(solutionTSP))
    for i in 1:length(solutionTSP)
        for j in 1:length(solutionTSP)
            tempsMinIJ[i,j] = Float64[]
            attenteIJ[i,j] = Float64[]
            operationMinIJ[i,j] = Int64[]
        end
    end
    for i in 1:length(solutionTSP)-1
        for j in i+1:length(solutionTSP)
            if j == i+1 # k = -1
                Tij, Aij = tempsOperationBi(instance, i, j, -1)
                operation = -1
                push!(tempsMinIJ[i, j], Tij)
                push!(attenteIJ[i, j], Aij)
                push!(operationMinIJ[i, j], operation)
            else
                # k > 0
                for k in (i+1):(j-1)
                    Tij, Aij = tempsOperationBi(instance, i, j, k)
                    operation = k
                    push!(tempsMinIJ[i, j], Tij)
                    push!(attenteIJ[i, j], Aij)
                    push!(operationMinIJ[i, j], operation)
                end
            end
        end
    end
    return tempsMinIJ, operationMinIJ, attenteIJ
end

########## Filtrage pendant (iterative dominance filtering) ##########

function AEPbiFiltrage(instance)
    tempsMinIJ, operationMinIJ, attenteIJ = calculMeilleuresOperationsBi(instance)
    tempsMinIJ, attenteIJ, operationMinIJ = filtrageOperations(tempsMinIJ, operationMinIJ, attenteIJ, instance.solution.ordreVisiteCamion) # Filtrage par dominance des opérations T(i,j)

    V = []
    A = []
    P = []
    for i in 1:length(instance.solution.ordreVisiteCamion)
        push!(V, Float64[])
        push!(A, Float64[])
        push!(P, [])
    end

    push!(V[1], 0.0)
    push!(A[1], 0.0)
    push!(P[1], [-1,-1,-1]) # A changer
    for i in 2:length(V)
        #println(" $(i-1) : $(length(V[i-1]))")
        for k in 1:(i-1)
            for y in 1:length(V[k])
                for z in 1:length(tempsMinIJ[k,i])
                    temps = V[k][y] + (tempsMinIJ[k,i])[z]
                    attente = A[k][y] + (attenteIJ[k,i])[z]
                    push!(V[i], temps)
                    push!(A[i], attente)
                    push!(P[i], [k,y, trunc(Int64, (operationMinIJ[k,i])[z]) ]) # On vient de V[k,m] et on a pris l'opération T(k,i,n)
                end
            end
        end
        Vtemp = []
        Atemp = []
        Ptemp = []
        for a in 1:length(V[i]) # On veut tester si le point i est dominé par un point j
            estDomine = false
            for b in 1:length(V[i])
                if b != a
                    if round(V[i][b], digits=5) == round(V[i][a], digits=5) && round(A[i][b], digits=5) == round(A[i][a], digits=5)
                        println(" Egalité V[i]")
                    end
                    if round(V[i][b], digits=5) <= round(V[i][a], digits=5) && round(A[i][b], digits=5) <= round(A[i][a], digits=5)
                        estDomine = true
                        break
                    end
                end
            end
            if estDomine == false
                push!(Vtemp, V[i][a])
                push!(Atemp, A[i][a])
                push!(Ptemp, P[i][a])
            end
        end
        V[i] = Vtemp
        A[i] = Atemp
        P[i] = Ptemp
    end
    #println(" V = $V")
    #println(" A = $A")
    #println(" P = $P")
    return V, A, P
end

########## Filtrage par dominance des T(i,j) ##########

function filtrageOperations(tempsMinIJ, operationMinIJ, attenteIJ, solutionTSP)
    tempsMinIJBis = Array{Vector{Float64},2}(undef, length(solutionTSP), length(solutionTSP))
    attenteIJBis = Array{Vector{Float64},2}(undef, length(solutionTSP), length(solutionTSP))
    operationMinIJBis = Array{Vector{Float64},2}(undef, length(solutionTSP), length(solutionTSP))
    for i in 1:length(solutionTSP)
        for j in 1:length(solutionTSP)
            tempsMinIJBis[i,j] = Float64[]
            attenteIJBis[i,j] = Float64[]
            operationMinIJBis[i,j] = Int64[]
        end
    end

    for i in 1:length(solutionTSP)
        for j in i+1:length(solutionTSP)
            for l1 in 1:length(tempsMinIJ[i,j]) # Pour toutes les opérations dans T(i,j)
                estDomine = false # Test du point (i,j,l1)
                for l2 in 1:length(tempsMinIJ[i,j]) # Pour toutes les opérations dans T(i,j)
                    if l1 != l2
                        if round(tempsMinIJ[i,j][l2], digits=5) == round(tempsMinIJ[i,j][l1], digits=5) && round(attenteIJ[i,j][l2], digits=5) == round(attenteIJ[i,j][l1], digits=5) # Cas d'égalité
                            println(" Egalité operation")
                        end
                        if round(tempsMinIJ[i,j][l2], digits=5) <= round(tempsMinIJ[i,j][l1], digits=5) && round(attenteIJ[i,j][l2], digits=5) <= round(attenteIJ[i,j][l1], digits=5)
                            estDomine = true
                            break
                        end
                    end
                end
                if estDomine == false
                    push!(tempsMinIJBis[i,j], tempsMinIJ[i,j][l1])
                    push!(attenteIJBis[i,j], attenteIJ[i,j][l1])
                    push!(operationMinIJBis[i,j], operationMinIJ[i,j][l1])
                end
            end
        end
    end
    return tempsMinIJBis, attenteIJBis, operationMinIJBis
end

########## Graham ##########

function orient(p,q,r)
    pr_x = r[1] - p[1]
    pr_y = r[2] - p[2]
    rq_x = q[1] - r[1]
    rq_y = q[2] - r[2]
    sign(pr_x * rq_y - pr_y * rq_x)
end

function convex_hull_graham_scan!(points)
    @assert length(points) >= 3 "Too few points"
    sort!(points)  # in lexical order

    function build_convex_hull(points)
        n = length(points)
        hull = points[1:2]
        sizehint!(hull, n)

        for i in 3:n
            while true
                if orient(hull[end - 1], hull[end], points[i]) > 0
                    break
                else
                    pop!(hull)
                    if length(hull) < 2
                        break
                    end
                end
            end
            push!(hull, points[i])
        end

        hull
    end

    upper_hull = build_convex_hull(points)
    pop!(upper_hull)
    lower_hull = build_convex_hull(@view points[end:-1:1])
    pop!(lower_hull)

    [upper_hull; lower_hull]
end

function convex_hull_length(hull)
    hull2 = circshift(hull, -1)
    map((x, y) -> hypot((x .- y)...), hull2, hull) |> sum
end

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver