text:
struct SortWeeks <: Constraint

    LeftGroup::Variable
    RightGroup::Variable

end

function SortWeeks(model, LeftGroup, RightGroup)
    constraint = SortWeeks(LeftGroup, RightGroup)

    push!(LeftGroup.linkedConstraint, constraint)
    push!(RightGroup.linkedConstraint, constraint)

    push!(model.constraints, constraint)
end

function filtrage!(weeks::SortWeeks)

    changeVariable = Vector{Variable}(undef,3)
    nbChange = 0

    changeLeft = Change(weeks.LeftGroup)
    changeRight = Change(weeks.RightGroup)

    #### filtering LeftGroup

    if !weeks.LeftGroup.isFixed
        # v : element max de RightGroup.upperBound
        v = maximum(weeks.RightGroup.upperBound)

        c = 0
        n_UB = copy(weeks.LeftGroup.upperBound)
        changed = false

        for e in 1:length(weeks.LeftGroup.upperBound)
            if weeks.LeftGroup.upperBound[e] >= v
                deleteat!(n_UB, e-c)
                c += 1

                push!(changeLeft.removed, weeks.LeftGroup.upperBound[e])
                changed = true
            end
        end

        if changed
            nbChange += 1
            changeVariable[nbChange] = weeks.LeftGroup
        end

        weeks.LeftGroup.upperBound = n_UB

        cardSup = min(weeks.LeftGroup.cardinalSup, length(weeks.LeftGroup.upperBound))
        changeLeft.cardRemoved = cardSup - weeks.LeftGroup.cardinalSup

        weeks.LeftGroup.cardinalSup += changeLeft.cardRemoved
    end

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

    #### filtering RightGroup
    if length(weeks.LeftGroup.lowerBound) >0 && length(weeks.RightGroup.lowerBound) > 0
        v = maximum(weeks.LeftGroup.lowerBound)
        w = maximum(weeks.RightGroup.lowerBound)

        if !weeks.RightGroup.isFixed && weeks.RightGroup.cardinalSup - length(weeks.RightGroup.lowerBound) == 1 && v > w

            c = 0
            n_UB = copy(weeks.RightGroup.upperBound)
            changed = false
            for e in 1:length(weeks.RightGroup.upperBound)
                if weeks.RightGroup.upperBound[e] < v && !(weeks.RightGroup.upperBound[e] in weeks.RightGroup.lowerBound)
                    deleteat!(n_UB, e-c)
                    c += 1

                    push!(changeRight.removed, weeks.RightGroup.upperBound[e])
                    changed = true
                end
            end

            if changed
                nbChange += 1
                changeVariable[nbChange] = weeks.RightGroup
            end

            weeks.RightGroup.upperBound = n_UB

            cardSup = min(weeks.RightGroup.cardinalSup, length(weeks.RightGroup.upperBound))
            changeRight.cardRemoved = cardSup - weeks.RightGroup.cardinalSup

            weeks.RightGroup.cardinalSup += changeRight.cardRemoved
        end
    end

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

    stop = weeks.LeftGroup.cardinalInf > weeks.LeftGroup.cardinalSup
	stop |= weeks.RightGroup.cardinalInf > weeks.RightGroup.cardinalSup

    return changeVariable[1:nbChange], (changeLeft, changeRight), stop
end

function quoi()
    v1 = Variable("v1", [1,12], [1,3,4,6,10,12,14], 4, 4, Set([]), false)
    v2 = Variable("v2", [1,5,10], [1,5,9,10,11,13,16], 4, 4, Set([]), false)
    println(v1.upperBound)
    println(v2.upperBound)
    model = ModelTest()
    sg = SortWeeks(v1,v2)
    println(sg)
    filtrage!(sg)
    println(v1.upperBound)
    println(v2.upperBound)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Serialization

import Base.MPFR
@testset "constructors" begin
    setprecision(53) do
        x = BigFloat()
        x = BigFloat(12)
    end
    x = BigFloat(12)
    @test x == BigFloat(x) == BigFloat(0xc) == BigFloat(12.) ==
          BigFloat(BigInt(12)) == BigFloat(BigFloat(12)) == parse(BigFloat,"12") ==
          parse(BigFloat,"12 ") == parse(BigFloat," 12") == parse(BigFloat," 12 ") ==
          BigFloat(Float32(12.)) == BigFloat(12//1) == BigFloat(SubString("12"))

    @test typeof(BigFloat(typemax(Int8))) == BigFloat
    @test typeof(BigFloat(typemax(Int16))) == BigFloat
    @test typeof(BigFloat(typemax(Int32))) == BigFloat
    @test typeof(BigFloat(typemax(Int64))) == BigFloat
    @test typeof(BigFloat(typemax(Int128))) == BigFloat

    @test typeof(BigFloat(true)) == BigFloat
    @test typeof(BigFloat(typemax(UInt8))) == BigFloat
    @test typeof(BigFloat(typemax(UInt16))) == BigFloat
    @test typeof(BigFloat(typemax(UInt32))) == BigFloat
    @test typeof(BigFloat(typemax(UInt64))) == BigFloat
    @test typeof(BigFloat(typemax(UInt128))) == BigFloat

    @test typeof(BigFloat(floatmax(Float32))) == BigFloat
    @test typeof(BigFloat(floatmax(Float64))) == BigFloat

    @test typeof(BigFloat(BigInt(1))) == BigFloat
    @test typeof(BigFloat(BigFloat(1))) == BigFloat

    @test typeof(BigFloat(1//1)) == BigFloat
    @test typeof(BigFloat(one(Rational{BigInt}))) == BigFloat
end
@testset "basic arithmetic" begin
    tol = 1e-12

    a = parse(BigFloat,"12.34567890121")
    b = parse(BigFloat,"12.34567890122")

    @test a+1e-11 ≈ b atol=tol
    @test !(b == a)
    @test b > a
    @test b >= a
    @test !(b < a)
    @test !(b <= a)

    c = parse(BigFloat,"24.69135780242")
    @test typeof(a * 2) == BigFloat
    @test a*2 ≈ c atol=tol
    @test (c-a) ≈ a atol=tol

    d = parse(BigFloat,"-24.69135780242")
    @test typeof(d) == BigFloat
    @test d+c ≈ 0 atol=tol

    @test (BigFloat(3)/BigFloat(2)) ≈ BigFloat(1.5) atol=tol
    @testset "+" begin
        x = BigFloat(12)
        y = BigFloat(30)
        @test x + y == BigFloat(42)
        @test x + typemax(UInt128) == x + BigInt(typemax(UInt128))
        @test x + typemax(Int128) == x + BigInt(typemax(Int128))
        f = parse(BigFloat,"1234567890.123")
        g = parse(BigFloat,"1234567891.123")

        tol = 1e-3

        @test f+Int8(1) ≈ g atol=tol
        @test f+Int16(1) ≈ g atol=tol
        @test f+Int32(1) ≈ g atol=tol
        @test f+Int64(1) ≈ g atol=tol
        @test f+Int128(1) ≈ g atol=tol

        @test f+true ≈ g atol=tol
        @test f+UInt8(1) ≈ g atol=tol
        @test f+UInt16(1) ≈ g atol=tol
        @test f+UInt32(1) ≈ g atol=tol
        @test f+UInt64(1) ≈ g atol=tol
        @test f+UInt128(1) ≈ g atol=tol

        @test f+BigInt(1) ≈ g atol=tol

        @test f+1f0 ≈ g atol=tol
        @test f+1e0 ≈ g atol=tol

        @test f+BigFloat(1) ≈ g atol=tol

        @test f+(1//1) ≈ g atol=tol

        @test f+one(Rational{BigInt}) ≈ g atol=tol

        a = parse(BigFloat,"123456789012345678901234567890")
        b = parse(BigFloat,"123456789012345678901234567891")
        @testset "Signed addition" begin
            @test a+Int8(1) == b
            @test a+Int16(1) == b
            @test a+Int32(1) == b
            @test a+Int64(1) == b
            @test Int8(1)+ a == b
            @test Int16(1)+a == b
            @test Int32(1)+a == b
            @test Int64(1)+a == b
            @test b+Int8(-1) == a
            @test b+Int16(-1) == a
            @test b+Int32(-1) == a
            @test b+Int64(-1) == a
            @test Int8(-1)+ b == a
            @test Int16(-1)+b == a
            @test Int32(-1)+b == a
            @test Int64(-1)+b == a
        end
        @testset "Unsigned addition" begin
            @test a+true == b
            @test a+UInt8(1) == b
            @test a+UInt16(1) == b
            @test a+UInt32(1) == b
            @test a+UInt64(1) == b
            @test true+a == b
            @test UInt8(1)+ a == b
            @test UInt16(1)+a == b
            @test UInt32(1)+a == b
            @test UInt64(1)+a == b
        end
        @testset "Float64 addition" begin
            @test a + 1.0f0 == b
            @test 1.0f0 + a == b
            @test a + 1.0 == b
            @test 1.0 + a == b
        end
        @testset "BigInt addition" begin
            @test a + BigInt(1) == b
            @test BigInt(1) + a == b
        end
    end
    @testset "-" begin
        x = BigFloat(12)
        y = BigFloat(-30)
        @test x - y == BigFloat(42)
        @test x - typemax(UInt128) == x - BigInt(typemax(UInt128))
        @test x - typemax(Int128) == x - BigInt(typemax(Int128))
        a = parse(BigFloat,"123456789012345678901234567890")
        b = parse(BigFloat,"123456789012345678901234567891")

        @testset "Signed subtraction" begin
            @test b-Int8(1) == a
            @test b-Int16(1) == a
            @test b-Int32(1) == a
            @test b-Int64(1) == a
            @test Int8(1)- b == -a
            @test Int16(1)-b == -a
            @test Int32(1)-b == -a
            @test Int64(1)-b == -a
            @test a-Int8(-1) == b
            @test a-Int16(-1) == b
            @test a-Int32(-1) == b
            @test a-Int64(-1) == b
            @test Int8(-1)- a == -b
            @test Int16(-1)-a == -b
            @test Int32(-1)-a == -b
            @test Int64(-1)-a == -b
        end
        @testset "Unsigned subtraction" begin
            @test b-true == a
            @test b-UInt8(1) == a
            @test b-UInt16(1) == a
            @test b-UInt32(1) == a
            @test b-UInt64(1) == a
            @test true-b == -a
            @test UInt8(1)- b == -a
            @test UInt16(1)-b == -a
            @test UInt32(1)-b == -a
            @test UInt64(1)-b == -a
        end
        @testset "Float64 subtraction" begin
            @test b - 1.0f0 == a
            @test 1.0f0 - b == -a
            @test b - 1.0 == a
            @test 1.0 - b == -a
        end
        @testset "BigInt subtraction" begin
            @test b - BigInt(1) == a
            @test BigInt(1) - b == -a
        end
    end
    @testset "*" begin
        x = BigFloat(6)
        y = BigFloat(9)
        @test x * y != BigFloat(42)
        @test x * y == BigFloat(54)
        @test x * typemax(UInt128) == x * BigInt(typemax(UInt128))
        @test x * typemax(Int128) == x * BigInt(typemax(Int128))

        a = parse(BigFloat,"123456789012345678901234567890")
        b = parse(BigFloat,"123456789012345678901234567891")
        @testset "Signed multiplication" begin
            @test a*Int8(1) == a
            @test a*Int16(1) == a
            @test a*Int32(1) == a
            @test a*Int64(1) == a
            @test Int8(1)* a == a
            @test Int16(1)*a == a
            @test Int32(1)*a == a
            @test Int64(1)*a == a
            @test a*Int8(-1) == -a
            @test a*Int16(-1) == -a
            @test a*Int32(-1) == -a
            @test a*Int64(-1) == -a
            @test Int8(-1)* a == -a
            @test Int16(-1)*a == -a
            @test Int32(-1)*a == -a
            @test Int64(-1)*a == -a
        end
        @testset "Unsigned multiplication" begin
            @test a*true == a
            @test a*UInt8(1) == a
            @test a*UInt16(1) == a
            @test a*UInt32(1) == a
            @test a*UInt64(1) == a
            @test true*a == a
            @test UInt8(1)* a == a
            @test UInt16(1)*a == a
            @test UInt32(1)*a == a
            @test UInt64(1)*a == a
        end
        @testset "Float64 multiplication" begin
            @test a * 1.0f0 == a
            @test 1.0f0 * a == a
            @test a * 1.0 == a
            @test 1.0 * a == a
        end
        @testset "BigInt multiplication" begin
            @test a * BigInt(1) == a
            @test BigInt(1) * a == a
        end
    end
    @testset "/" begin
        x = BigFloat(9)
        y = BigFloat(6)
        @test x / y == BigFloat(9/6)
        @test x / typemax(UInt128) == x / BigInt(typemax(UInt128))
        @test x / typemax(Int128) == x / BigInt(typemax(Int128))

        a = parse(BigFloat,"123456789012345678901234567890")
        b = parse(BigFloat,"123456789012345678901234567891")
        c = parse(BigInt,"61728394506172839450617283945")
        # d = 2^200
        d = parse(BigFloat,"1606938044258990275541962092341162602522202993782792835301376")
        f = parse(BigFloat,"6.223015277861141707144064053780124240590252168721167133101116614789698834035383e-61")

        @testset "Signed division" begin
            @test a/Int8(2) == c
            @test a/Int16(2) == c
            @test a/Int32(2) == c
            @test a/Int64(2) == c
            @test Int8(1)/ d == f
            @test Int16(1)/d == f
            @test Int32(1)/d == f
            @test Int64(1)/d == f
            @test a/Int8(-2) == -c
            @test a/Int16(-2) == -c
            @test a/Int32(-2) == -c
            @test a/Int64(-2) == -c
            @test Int8(-1)/ d == -f
            @test Int16(-1)/d == -f
            @test Int32(-1)/d == -f
            @test Int64(-1)/d == -f
        end
        @testset "Unsigned division" begin
            @test a/true == a
            @test a/UInt8(2) == c
            @test a/UInt16(2) == c
            @test a/UInt32(2) == c
            @test a/UInt64(2) == c
            @test true/d == f
            @test UInt8(1)/ d == f
            @test UInt16(1)/d == f
            @test UInt32(1)/d == f
            @test UInt64(1)/d == f
        end
        @testset "Float64 division" begin
            @test a / 2.0f0 == c
            @test 1.0f0 / d == f
            @test a / 2.0 == c
            @test 1.0 / d == f
        end
        # BigInt division
        @test a / BigInt(2) == c
        # inv
        @test inv(x) == one(x)/x == 1/x == x^-1 == Clong(1)/x
    end
    #^
    x = BigFloat(12)
    y = BigFloat(4)
    @test x^y == BigFloat(20736)

    @test big(2.0)^big(3) == 8
    for T in [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt]
        @test T(2)^big(3.0) == 8
        @test big(2.0)^T(3) == 8
    end
end
@testset "iterated arithmetic" begin
    a = BigFloat(12.25)
    b = BigFloat(23.125)
    c = BigFloat(-7)
    d = BigFloat(-12.75)
    f = BigFloat(2.0625)
    g = BigFloat(0.03125)
    @test +(a, b) == BigFloat(35.375)
    @test +(a, b, c) == BigFloat(28.375)
    @test +(a, b, c, d) == BigFloat(15.625)
    @test +(a, b, c, d, f) == BigFloat(17.6875)
    @test +(a, b, c, d, f, g) == BigFloat(17.71875)

    @test *(a, b) == parse(BigFloat,"2.8328125e+02")
    @test *(a, b, c) == parse(BigFloat,"-1.98296875e+03")
    @test *(a, b, c, d) == parse(BigFloat,"2.52828515625e+04")
    @test *(a, b, c, d, f) == parse(BigFloat,"5.214588134765625e+04")
    @test *(a, b, c, d, f, g) == parse(BigFloat,"1.6295587921142578125e+03")
end
@testset "< / > / <= / >=" begin
    x = BigFloat(12)
    y = BigFloat(42)
    z = BigFloat(30)
    @test y > x
    @test y >= x
    @test y > z
    @test y >= z
    @test x < y
    @test x <= y
    @test z < y
    @test z <= y
    @test y - x >= z
    @test y - x <= z
    @test !(x >= z)
    @test !(y <= z)
end
@testset "rounding modes" begin
    setprecision(4) do
        # default mode is round to nearest
        down, up =  setrounding(BigFloat,RoundNearest) do
            parse(BigFloat,"0.0938"), parse(BigFloat,"0.102")
        end
        setrounding(BigFloat,RoundDown) do
            @test BigFloat(0.1) == down
            @test BigFloat(0.1) != up
        end
        setrounding(BigFloat,RoundUp) do
            @test BigFloat(0.1) != down
            @test BigFloat(0.1) == up
        end
    end
end

@testset "copysign" begin
    x = BigFloat(1)
    y = BigFloat(-1)
    @test copysign(x, y) == y
    @test copysign(y, x) == x
    @test copysign(1.0, BigFloat(NaN)) == 1.0
end
@testset "isfinite / isinf / isnan" begin
    x = BigFloat(Inf)
    y = BigFloat(1)
    @test isinf(x) == true
    @test isinf(y) == false
    @test isfinite(x) == false
    @test isinf(x) == true

    x = BigFloat(NaN)
    y = BigFloat(1)
    @test isnan(x) == true
    @test isnan(y) == false
end
@testset "convert to BigFloat" begin
    @test convert(BigFloat, 1//2) == parse(BigFloat,"0.5")
    @test typeof(convert(BigFloat, 1//2)) == BigFloat
    @test convert(BigFloat, 0.5) == parse(BigFloat,"0.5")
    @test typeof(convert(BigFloat, 0.5)) == BigFloat
    @test convert(BigFloat, 40) == parse(BigFloat,"40")
    @test typeof(convert(BigFloat, 40)) == BigFloat
    @test convert(BigFloat, Float32(0.5)) == parse(BigFloat,"0.5")
    @test typeof(convert(BigFloat, Float32(0.5))) == BigFloat
    @test convert(BigFloat, parse(BigInt,"9223372036854775808")) == parse(BigFloat,"9223372036854775808")
    @test typeof(convert(BigFloat, parse(BigInt,"9223372036854775808"))) == BigFloat
    @test convert(AbstractFloat, parse(BigInt,"9223372036854775808")) == parse(BigFloat,"9223372036854775808")
    @test typeof(convert(AbstractFloat, parse(BigInt,"9223372036854775808"))) == BigFloat

    @test signbit(BigFloat(NaN)) == 0
    @test signbit(BigFloat(-NaN)) == 1
end
@testset "convert from BigFloat" begin
    @test convert(Float64, BigFloat(0.5)) == 0.5
    @test convert(Float32, BigFloat(0.5)) == Float32(0.5)
    @test convert(Float16, BigFloat(0.5)) == Float16(0.5)
    @test convert(Bool, BigFloat(0.0)) == false
    @test convert(Bool, BigFloat(1.0)) == true
    @test_throws InexactError convert(Bool, BigFloat(0.1))

    @test signbit(Float64(BigFloat(NaN))) == 0
    @test signbit(Float64(-BigFloat(NaN))) == 1
end
@testset "exponent, frexp, significand" begin
    x = BigFloat(0)
    @test_throws DomainError exponent(x)
    x = BigFloat(Inf)
    @test_throws DomainError exponent(x)
    x = BigFloat(15.674)
    @test exponent(x) == exponent(15.674)

    for i in [big(0.2), big(1.2), big(1220.0), big(23414.123)]
        mantissa, ex = frexp(i)
        @test i == mantissa * 2. ^ ex
    end

    for i in [big(0.2), big(1.2), big(1220.0), big(23414.123)]
        @test i == significand(i) * 2. ^ exponent(i)
    end
end
@testset "nextfloat/prevfloat" begin
    #should be immutable
    x = 12.
    y = BigFloat(x)
    @test x == y
    nextfloat(y)
    @test x == y
    prevfloat(y)
    @test x == y
    setprecision(53) do
        x = BigFloat(12.12)
        @test BigFloat(nextfloat(12.12)) == nextfloat(x)
        @test BigFloat(prevfloat(12.12)) == prevfloat(x)
    end
    @test isnan(nextfloat(BigFloat(NaN)))
    @test isnan(prevfloat(BigFloat(NaN)))
end
# sqrt DomainError
@test_throws DomainError sqrt(BigFloat(-1))

@testset "precision" begin
    old_precision = precision(BigFloat)
    x = BigFloat(0)
    @test precision(x) == old_precision
    setprecision(256)
    x = BigFloat(0)
    @test precision(x) == 256
    setprecision(old_precision)
    z = setprecision(240) do
        z = x + 20
        return z
    end
    @test float(z) == 20.
    @test precision(z) == 240
    x = BigFloat(12)
    @test precision(x) == old_precision
    @test_throws DomainError setprecision(1)
    # issue 15659
    @test (setprecision(53) do; big(1/3); end) < 1//3
end
@testset "isinteger" begin
    @test !isinteger(BigFloat(1.2))
    @test isinteger(BigFloat(12))
    @test isinteger(zero(BigFloat))
    @test isinteger(-zero(BigFloat))
    @test !isinteger(nextfloat(zero(BigFloat)))
    @test !isinteger(prevfloat(zero(BigFloat)))
    @test isinteger(maxintfloat(BigFloat))
    @test isinteger(-maxintfloat(BigFloat))
    @test !isinteger(BigFloat(Inf))
    @test !isinteger(-BigFloat(Inf))
    @test !isinteger(BigFloat(NaN))
end

@testset "comparisons" begin
    x = BigFloat(1)
    y = BigFloat(-1)
    z = BigFloat(NaN)
    ipl = BigFloat(Inf)
    imi = BigFloat(-Inf)
    @test x > y
    @test x >= y
    @test x >= x
    @test y < x
    @test y <= x
    @test y <= y
    @test x < ipl
    @test x <= ipl
    @test x > imi
    @test x >= imi
    @test imi == imi
    @test ipl == ipl
    @test imi < ipl
    @test z != z
    @test !(z == z)
    @test !(z <= z)
    @test !(z < z)
    @test !(z >= z)
    @test !(z > z)

    @test !isequal(BigFloat(0.0),BigFloat(-0.0))
    @test isequal(z, BigFloat(NaN))

    # total ordering
    @test isless(big(-0.0), big(0.0))
    @test isless(big(1.0), big(NaN))
    @test isless(big(Inf), big(NaN))
    @test isless(big(Inf), -big(NaN))
    @test !isless(big(NaN), big(NaN))
    @test !isless(big(-NaN), big(NaN))
    @test !isless(-big(NaN), big(NaN))
    @test !isless(-big(NaN), big(1.0))
    @test !isless(-big(NaN), 1.0)

    # cmp
    @test cmp(big(-0.0), big(0.0)) == -1
    @test cmp(big(0.0), big(-0.0)) == 1
    @test cmp(big(1.0), big(NaN)) == -1
    @test cmp(big(NaN), big(NaN)) == 0
    @test cmp(big(NaN), big(1.0)) == 1
end
@testset "signbit" begin
    @test signbit(BigFloat(-1.0)) == 1
    @test signbit(BigFloat(1.0)) == 0
    @test signbit(BigFloat(-0.0)) == 1
end
@testset "modf" begin
    x = BigFloat(12)
    y = BigFloat(0.5)
    @test modf(x+y) == (y, x)
    x = BigFloat(NaN)
    @test map(isnan, modf(x)) == (true, true)
    x = BigFloat(Inf)
    y = modf(x)
    @test (isnan(y[1]), isinf(y[2])) == (true, true)
end
@testset "rem" begin
    setprecision(53) do
        x = BigFloat(2)
        y = BigFloat(1.67)
        @test rem(x,y) == rem(2, 1.67)
        y = BigFloat(NaN)
        @test isnan(rem(x,y))
        @test isnan(rem(y,x))
        y = BigFloat(Inf)
        @test rem(x,y) == x
        @test isnan(rem(y,x))
    end
end
@testset "min/max" begin
    x = BigFloat(4)
    y = BigFloat(2)
    @test max(x,y) == x
    @test min(x,y) == y
    y = BigFloat(NaN)
    @test isnan(max(x,y))
    @test isnan(min(x,y))
    @test isnan(max(y,y))
    @test isnan(min(y,y))
end
@testset "sum" begin
    x = BigFloat(1)
    y = BigFloat(2)
    z = BigFloat(3)
    w = BigFloat(4)
    @test sum([x,y,z,w]) == BigFloat(10)
    big_array = fill(BigFloat(1), 100)
    @test sum(big_array) == BigFloat(100)
    @test sum(BigFloat[]) == BigFloat(0)
end
@testset "promotion" begin
    # the array converts everyone to the DEFAULT_PRECISION!
    x = BigFloat(12)
    y = setprecision(60) do
        BigFloat(42)
    end
    @test [x,y] == [BigFloat(12), BigFloat(42)]
end
@testset "log / log2 / log10" begin
    setprecision(53) do
    x = BigFloat(42)
        @test log(x) == log(42)
        @test isinf(log(BigFloat(0)))
        @test_throws DomainError log(BigFloat(-1))
        @test log2(x) == log2(42)
        @test isinf(log2(BigFloat(0)))
        @test_throws DomainError log2(BigFloat(-1))
        @test log10(x) == log10(42)
        @test isinf(log10(BigFloat(0)))
        @test_throws DomainError log10(BigFloat(-1))
    end
end
@testset "exp / exp2 / exp10" begin
    setprecision(53) do
        x = BigFloat(10)
        @test exp(x) == exp(10)
        @test exp2(x) == 1024
        @test exp10(x) == 10000000000
    end
end
@testset "convert to integer types" begin
    x = BigFloat(12.1)
    y = BigFloat(42)
    @test_throws InexactError convert(Int32, x)
    @test_throws InexactError convert(Int64, x)
    @test_throws InexactError convert(BigInt, x)
    @test_throws InexactError convert(UInt32, x)
    @test_throws InexactError convert(UInt32, x)
    @test convert(Int32, y) == 42
    @test convert(Int64, y) == 42
    @test convert(BigInt, y) == 42
    @test convert(UInt32, y) == 42
    @test convert(UInt32, y) == 42
end
@testset "round" begin
    x = BigFloat(42.42)
    y = setprecision(256) do
        parse(BigFloat,"9223372036854775809.2324")
    end
    z = parse(BigInt,"9223372036854775809")
    @test round(Integer,x) == 42
    @test round(Integer,y) == z
    @test typeof(round(UInt8, x)) == UInt8 && round(UInt8, x) == 0x2a
    @test typeof(round(UInt16, x)) == UInt16 && round(UInt16, x) == 0x2a
    @test typeof(round(UInt32, x)) == UInt32 && round(UInt32, x) == 0x2a
    @test typeof(round(UInt64, x)) == UInt64 && round(UInt64, x) == 0x2a
    @test typeof(round(Int64, x)) == Int64 && round(Int64, x) == 42
    @test typeof(round(Int, x)) == Int && round(Int, x) == 42
    @test typeof(round(UInt, x)) == UInt && round(UInt, x) == 0x2a
end
@testset "string representation" begin
    str = "1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012"
    setprecision(406) do
        @test string(nextfloat(BigFloat(1))) == str
    end
    setprecision(21) do
        @test string(parse(BigFloat, "0.1")) == "1.0000002e-01"
        @test string(parse(BigFloat, "-9.9")) == "-9.9000015"
    end
    setprecision(40) do
        @test string(parse(BigFloat, "0.1")) == "1.0000000000002e-01"
        @test string(parse(BigFloat, "-9.9")) == "-9.8999999999942"
    end
    setprecision(123) do
        @test string(parse(BigFloat, "0.1")) == "9.99999999999999999999999999999999999953e-02"
        @test string(parse(BigFloat, "-9.9")) == "-9.8999999999999999999999999999999999997"
    end
end
@testset "eps" begin
    x = eps(BigFloat)
    @test BigFloat(1) + x == BigFloat(1) + prevfloat(x)
    @test eps(BigFloat) == eps(BigFloat(1))
end
@testset "floatmin/floatmax" begin
    x = floatmin(BigFloat)
    @test x > 0
    @test prevfloat(x) == 0
    x = floatmax(BigFloat)
    @test !isinf(x)
    @test isinf(nextfloat(x))
end
@testset "factorial" begin
    setprecision(256) do
        x = BigFloat(42)
        @test factorial(x) == factorial(BigInt(42))
        x = BigFloat(10)
        @test factorial(x) == factorial(10)
        @test_throws DomainError factorial(BigFloat(-1))
        @test_throws DomainError factorial(BigFloat(331.3))
    end
end
@testset "trigonometric functions" begin
    setprecision(53) do
        for f in (:sin,:cos,:tan,:sec,:csc,:cot,:acos,:asin,:atan,
                :cosh,:sinh,:tanh,:sech,:csch,:coth,:asinh),
            j in (-1., -0.5, -0.25, .25, .5, 1.)
            @eval begin
                @test ($f)(BigFloat($j)) ≈ ($f)($j)
            end
        end
        for f in (:acos,:asin,:acosh,:atanh),
            j in (-2, -1.5)
            @eval begin
                @test_throws DomainError ($f)(BigFloat($j))
            end
        end
        for f in (:sin,:cos,:tan,:sec,:csc,:cot,:cosh,:sinh,:tanh,
                  :sech,:csch,:coth,:acosh,:asinh),
            j in (1., 1.5, 1.9)
            @eval begin
                @test ($f)(BigFloat($j)) ≈ ($f)($j)
            end
        end
        for j in (.25, .5)
            @test atanh(BigFloat(j)) ≈ atanh(j)
        end
    end

    # hypot
    @test hypot(BigFloat(3), BigFloat(4)) == 5

    # atan
    setprecision(53) do
        @test atan(12,2) == atan(BigFloat(12), BigFloat(2))
    end
end
@testset "ldexp" begin
    setprecision(53) do
        @test ldexp(BigFloat(24.5), 72) == ldexp(24.5, 72)
        @test ldexp(BigFloat(24.5), Int16(72)) == ldexp(24.5, 72)
        @test ldexp(BigFloat(24.5), -72) == ldexp(24.5, -72)
        @test ldexp(BigFloat(24.5), Int16(-72)) == ldexp(24.5, -72)
        @test ldexp(BigFloat(24.5), UInt(72)) == ldexp(24.5, 72)
        @test ldexp(BigFloat(24.5), 0x48) == ldexp(24.5, 72)
    end
end
@testset "ceil / floor / trunc" begin
    x = BigFloat(12.042)
    @test BigFloat(13) == ceil(x)
    x = parse(BigFloat,"28273.7312487489135135135")
    y = BigInt(28273)
    z = BigInt(28274)
    a = parse(BigFloat,"123456789012345678901234567890.2414")
    b = parse(BigInt,"123456789012345678901234567890")
    c = parse(BigInt,"123456789012345678901234567891")
    @test ceil(x) == z
    @test typeof(ceil(x)) == BigFloat
    @test floor(x) == y
    @test typeof(floor(x)) == BigFloat
    @test trunc(x) == y
    @test typeof(trunc(x)) == BigFloat

    @test ceil(Integer,x) == z
    @test typeof(ceil(Integer,x)) == BigInt
    @test floor(Integer,x) == y
    @test typeof(floor(Integer,x)) == BigInt
    @test trunc(Integer,x) == y
    @test typeof(trunc(Integer,x)) == BigInt

    @test ceil(Int64, x) == Int64(z)
    @test typeof(ceil(Int64, x)) == Int64
    @test floor(Int64, x) == Int64(y)
    @test typeof(floor(Int64, x)) == Int64
    @test trunc(Int64, x) == Int64(y)
    @test typeof(trunc(Int64, x)) == Int64

    @test ceil(Int32, x) == Int32(z)
    @test typeof(ceil(Int32, x)) == Int32
    @test floor(Int32, x) == Int32(y)
    @test typeof(floor(Int32, x)) == Int32
    @test trunc(Int32, x) == Int32(y)
    @test typeof(trunc(Int32, x)) == Int32

    @test ceil(Int16, x) == Int16(z)
    @test typeof(ceil(Int16, x)) == Int16
    @test floor(Int16, x) == Int16(y)
    @test typeof(floor(Int16, x)) == Int16
    @test trunc(Int16, x) == Int16(y)
    @test typeof(trunc(Int16, x)) == Int16

    #@test ceil(Int8, x) == Int8(z)
    #@test typeof(ceil(Int8, x)) == Int8
    #@test floor(Int8, x) == Int8(y)
    #@test typeof(floor(Int8, x)) == Int8
    #@test trunc(Int8, x) == Int8(y)
    #@test typeof(trunc(Int8, x)) == Int8

    @test ceil(UInt64, x) == UInt64(z)
    @test typeof(ceil(UInt64, x)) == UInt64
    @test floor(UInt64, x) == UInt64(y)
    @test typeof(floor(UInt64, x)) == UInt64
    @test trunc(UInt64, x) == UInt64(y)
    @test typeof(trunc(UInt64, x)) == UInt64

    @test ceil(UInt32, x) == UInt32(z)
    @test typeof(ceil(UInt32, x)) == UInt32
    @test floor(UInt32, x) == UInt32(y)
    @test typeof(floor(UInt32, x)) == UInt32
    @test trunc(UInt32, x) == UInt32(y)
    @test typeof(trunc(UInt32, x)) == UInt32

    @test ceil(UInt16, x) == UInt16(z)
    @test typeof(ceil(UInt16, x)) == UInt16
    @test floor(UInt16, x) == UInt16(y)
    @test typeof(floor(UInt16, x)) == UInt16
    @test trunc(UInt16, x) == UInt16(y)
    @test typeof(trunc(UInt16, x)) == UInt16

    #@test ceil(UInt8, x) == UInt8(z)
    #@test typeof(ceil(UInt8, x)) == UInt8
    #@test floor(UInt8, x) == UInt8(y)
    #@test typeof(floor(UInt8, x)) == UInt8
    #@test trunc(UInt8, x) == UInt8(y)
    #@test typeof(trunc(UInt8, x)) == UInt8

    @test ceil(Integer,a) == c
    @test typeof(ceil(Integer,a)) == BigInt
    @test floor(Integer,a) == b
    @test typeof(floor(Integer,a)) == BigInt
    @test trunc(Integer,a) == b
    @test typeof(trunc(Integer,a)) == BigInt

    @test ceil(Int128,a) == c
    @test typeof(ceil(Int128,a)) == Int128
    @test floor(Int128,a) == b
    @test typeof(floor(Int128,a)) == Int128
    @test trunc(Int128,a) == b
    @test typeof(trunc(Int128,a)) == Int128

    @test ceil(UInt128,a) == c
    @test typeof(ceil(UInt128,a)) == UInt128
    @test floor(UInt128,a) == b
    @test typeof(floor(UInt128,a)) == UInt128
    @test trunc(UInt128,a) == b
    @test typeof(trunc(UInt128,a)) == UInt128
end
@testset "div" begin
    @test div(big"1.0",big"0.1") == 9
    @test div(1,big"0.1") == 9
    @test div(1.0,big"0.1") == 9
    @test div(big"1.0",0.1) == 9
    @test div(big"1",big"0.1") == 9
    @test div(big"1",0.1) == 9
end
@testset "issue #5963" begin
    @test typemax(Int128) == convert(BigFloat, typemax(Int128))
    @test typemax(Int128)  * big(1.0) == convert(BigFloat, typemax(Int128))
    @test typemax(UInt64)  * big(1.0) == big(typemax(UInt64))
    @test typemax(UInt128) * big(1.0) == big(typemax(UInt128))
end
@testset "issue #3399" begin
    i1 = BigInt(10)^Int32(1000)
    i2 = parse(BigInt,"10000000000000000000000000000000000000000000000000000000000000000000000000000040000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000")
    f = BigFloat(10)^Int32(1000)
    @test i1 != i2
    @test i1 != f
    @test i2 != f

    @test f > i1
    @test f > i2

    i3 = trunc(Integer,f)
    @test i3 == f
    @test i3+1 > f
    @test i3+1 >= f
end
# issue #8318
@test convert(Int64,big(500_000_000_000_000.)) == 500_000_000_000_000

@testset "issue #9816" begin
    # check exponent range is set to max possible
    @test MPFR.get_emin() == MPFR.get_emin_min()
    @test MPFR.get_emax() == MPFR.get_emax_max()
end
# issue #10994: handle embedded NUL chars for string parsing
@test_throws ArgumentError parse(BigFloat, "1\0")

@testset "serialization (issue #12386)" begin
    b = PipeBuffer()
    let x = setprecision(53) do
            return 2.1 * big(pi)
        end
        serialize(b, x)
        @test deserialize(b) == x
    end
    let x = BigFloat(Inf, 46)
        serialize(b, x)
        @test deserialize(b) == x == BigFloat(Inf, 2)
    end
end
@test isnan(sqrt(BigFloat(NaN)))

@testset "PR 17217" begin
    # PR 17217 -- BigFloat constructors with given precision and rounding mode
    # test constructors and `big` with additional precision and rounding mode:
    for prec in (10, 100, 1000)
        for val in ("3.1", pi, "-1.3", 3.1)
            let a = BigFloat(val),
                b = BigFloat(val, prec),
                c = BigFloat(val, RoundUp),
                d = BigFloat(val, prec, RoundDown),
                e = BigFloat(val, prec, RoundUp)

                @test precision(a) == precision(BigFloat)
                @test precision(b) == prec
                @test precision(c) == precision(BigFloat)
                @test precision(d) == prec
                @test precision(e) == prec
                (val != 3.1) && @test e > d     # rounding has no effect when constructing from Float64
            end
        end
    end
end

# issue #22758
if MPFR.version() > v"3.1.5" || "r11590" in MPFR.patches()
    setprecision(2_000_000) do
        @test abs(sin(big(pi)/6) - 0.5) < ldexp(big(1.0),-1_999_000)
    end
end

@testset "show BigFloat" begin
    function test_show_bigfloat(x::BigFloat; contains_e::Bool=true,
        ends::String="",
        starts::String="")
        sx = sprint(show, x)
        scx = sprint(show, x, context=:compact => true)
        strx = string(x)
        @test sx == strx
        @test length(scx) < 20
        @test length(scx) <= length(sx)
        @test x == parse(BigFloat, sx)
        @test ≈(x, parse(BigFloat, scx), rtol=1e-4)
        for s in (sx, scx)
            @test occursin('e', s) == contains_e
            @test startswith(s, starts)
            @test endswith(s, ends)
        end
    end

    test_show_bigfloat(big"1.23456789", contains_e=false, starts="1.23")
    test_show_bigfloat(big"-1.23456789", contains_e=false, starts="-1.23")
    test_show_bigfloat(big"2.3457645687563543266576889678956787e10000", starts="2.345", ends="e+10000")
    test_show_bigfloat(big"-2.3457645687563543266576889678956787e-10000", starts="-2.345", ends="e-10000")

    for to_string in [string,
                      x->sprint(show, x),
                      x->sprint(show, x, context=:compact => true)]
        @test to_string(big"0.0") == "0.0"
        @test to_string(big"-0.0") == "-0.0"
        @test to_string(big"1.0") == "1.0"
        @test to_string(big"-1.0") == "-1.0"
    end
end

@testset "big(::Type)" begin
    for x in (2f0, pi, 7.8, big(ℯ))
        @test big(typeof(x)) == typeof(big(x))
        @test big(typeof(complex(x, x))) == typeof(big(complex(x, x)))
    end
end

ext: jl

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
"""
    collect_top_level_nodes!
To build an expression DAG, we need to collect all of the symbols in advance.
"""
function collect_top_level_nodes! end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLCursor, options)
    dumpobj(cursor)
    file, line, col = get_file_line_column(cursor)
    error("No support for code at $file:$line:$col, please file an issue to Clang.jl.")
    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLFunctionDecl, options)
    func_type = getCursorType(cursor)

    if kind(func_type) == CXType_FunctionNoProto
        ty = FunctionNoProto()
    elseif isVariadic(func_type)
        ty = FunctionVariadic()
    elseif kind(func_type) == CXType_FunctionProto
        ty = FunctionProto()
    else
        ty = FunctionDefault()
    end

    id = Symbol(spelling(cursor))

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLTypedefDecl, options)
    lhs_type = getTypedefDeclUnderlyingType(cursor)

    if has_elaborated_reference(lhs_type)
        ty = TypedefElaborated()
    elseif has_function_reference(lhs_type)
        ty = TypedefFunction()
    else
        ty = TypedefDefault()
    end

    id = Symbol(spelling(cursor))

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLMacroDefinition, options)
    if isMacroBuiltin(cursor)
        ty = MacroBuiltIn()
    elseif isMacroFunctionLike(cursor)
        ty = MacroFunctionLike()
    elseif is_macro_definition_only(cursor)
        ty = MacroDefinitionOnly()
    elseif is_macro_unsupported(cursor)
        ty = MacroUnsupported()
    else
        ty = MacroDefault()
    end

    id = Symbol(spelling(cursor))

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLStructDecl, options)
    use_deterministic_sym = get(options, "use_deterministic_symbol", false)

    str = spelling(cursor)

    if isempty(str)
        @assert isCursorDefinition(cursor)
        ty = StructAnonymous()
        id = use_deterministic_sym ? gensym_deterministic("Ctag") : gensym("Ctag")
    elseif is_forward_declaration(cursor)
        ty = StructForwardDecl()
        id = Symbol(str)
    elseif isCursorDefinition(cursor)
        ty = StructDefinition()
        id = Symbol(str)
    else
        ty = StructDefault()
        id = Symbol(str)
    end

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLUnionDecl, options)
    use_deterministic_sym = get(options, "use_deterministic_symbol", false)

    str = spelling(cursor)

    if isempty(str)
        @assert isCursorDefinition(cursor)
        ty = UnionAnonymous()
        id = use_deterministic_sym ? gensym_deterministic("Ctag") : gensym("Ctag")
    elseif is_forward_declaration(cursor)
        ty = UnionForwardDecl()
        id = Symbol(str)
    elseif isCursorDefinition(cursor)
        ty = UnionDefinition()
        id = Symbol(str)
    else
        ty = UnionDefault()
        id = Symbol(str)
    end

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

function collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLEnumDecl, options)
    use_deterministic_sym = get(options, "use_deterministic_symbol", false)

    str = spelling(cursor)

    if isempty(str)
        @assert isCursorDefinition(cursor)
        ty = EnumAnonymous()
        id = use_deterministic_sym ? gensym_deterministic("Ctag") : gensym("Ctag")
    elseif is_forward_declaration(cursor)
        ty = EnumForwardDecl()
        id = Symbol(str)
    elseif isCursorDefinition(cursor)
        ty = EnumDefinition()
        id = Symbol(str)
    else
        ty = EnumDefault()
        id = Symbol(str)
    end

    push!(nodes, ExprNode(id, ty, cursor, Expr[], Int[]))

    return nodes
end

# skip macro expansion since the expanded info is already embedded in the AST
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLMacroInstantiation, options) = nodes
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLMacroExpansion, options) = nodes

# skip variable definition for now
# isn't it insane to define a variable in an interface header file? what's the use case?
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLVarDecl, options) = nodes

# skip `#include <...>`
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLInclusionDirective, options) = nodes
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLLastPreprocessing, options) = nodes  # FIXME: fix cltype.jl

# skip unexposed decl
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLUnexposedDecl, options) = nodes
collect_top_level_nodes!(nodes::Vector{ExprNode}, cursor::CLFirstDecl, options) = nodes  # FIXME: fix cltype.jl


text:
"""
    BetaPlane{T} <: AbstractRotation

A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`)
that accounts for the variation of the locally vertical component of the rotation
vector with latitude.
"""
struct BetaPlane{T} <: AbstractRotation
    f₀ :: T
     β :: T
end

"""
    BetaPlane([T=Float64;] f₀=nothing, β=nothing,
                           rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

The user may specify both `f₀` and `β`, or the three parameters `rotation_rate`,
`latitude`, and `radius` that specify the rotation rate and radius of a planet, and
the central latitude (where y = 0) at which the `β`-plane approximation is to be made.

By default, the `rotation_rate` and planet `radius` is assumed to be Earth's.
"""
function BetaPlane(T=Float64; f₀=nothing, β=nothing,
                              rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

    use_f_and_β = !isnothing(f₀) && !isnothing(β)
    use_planet_parameters = !isnothing(latitude)

    if !xor(use_f_and_β, use_planet_parameters)
        throw(ArgumentError("Either both keywords f₀ and β must be specified, " *
                            "*or* all of rotation_rate, latitude, and radius."))
    end

    if use_planet_parameters
        f₀ = 2rotation_rate * sind(latitude)
         β = 2rotation_rate * cosd(latitude) / radius
     end

    return BetaPlane{T}(f₀, β)
end

@inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) =
    @inbounds - (coriolis.f₀ + coriolis.β * grid.yC[j]) * ℑxyᶠᶜᵃ(i, j, k, grid, U[2])

@inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) =
    @inbounds   (coriolis.f₀ + coriolis.β * grid.yF[j]) * ℑxyᶜᶠᵃ(i, j, k, grid, U[1])

@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::BetaPlane, U) where FT = zero(FT)

Base.show(io::IO, β_plane::BetaPlane{FT}) where FT =
    print(io, "BetaPlane{$FT}: ", @sprintf("f₀ = %.2e, β = %.2e", β_plane.f₀, β_plane.β))


text:
# Lewis, E. E., et al. "Benchmark specification for Deterministic 2-D/3-D MOX fuel
# assembly transport calculations without spatial homogenization (C5G7 MOX)."
# NEA/NSC 280 (2001): 2001.
using MOCNeutronTransport

filename = "c5g7.step"
add_timestamps_to_logger()

gmsh.initialize()
gmsh.option.set_number("General.NumThreads", 0) # 0 uses system default, i.e. OMP_NUM_THREADS)
gmsh.option.set_number("Geometry.OCCParallel", 1) # use parallel OCC boolean operations
gmsh.option.set_number("General.Verbosity", 2) # 1: +errors, 2: +warnings, 3: +direct, 4: +information, 5: +status, 99: +debug

materials = import_model(filename, names = true)

lat_div = Vec{4,Float64}([i*21.42 for i = 0:3])
lattice_grid = RectilinearGrid(lat_div, lat_div)
module_grid = lattice_grid # assembly modular ray tracing
coarse_div = Vec{52,Float64}([1.26*i for i ∈ 0:17*3])
coarse_grid = RectilinearGrid(coarse_div, coarse_div)
mpact_grid = MPACTGridHierarchy(lattice_grid, module_grid, coarse_grid)

# See http://juliagraphics.github.io/Colors.jl/stable/namedcolors/ for named colors
# or use an RGB tuple
push!(materials, Material(name="Moderator", color="sienna", mesh_size=1.0))
overlay_mpact_grid_hierarchy(mpact_grid, materials)

Σₜᵢ_FC  = [1.90730E-01, 4.56520E-01, 6.40700E-01, 6.49840E-01, 6.70630E-01, 8.75060E-01, 1.43450E+00]
Σₜᵢ_GT  = [1.90730E-01, 4.56520E-01, 6.40670E-01, 6.49670E-01, 6.70580E-01, 8.75050E-01, 1.43450E+00]
Σₜᵢ_M43 = [2.11920E-01, 3.55810E-01, 4.88900E-01, 5.71940E-01, 4.32390E-01, 6.84950E-01, 6.88910E-01]
Σₜᵢ_M70 = [2.14540E-01, 3.59350E-01, 4.98910E-01, 5.96220E-01, 4.80350E-01, 8.39360E-01, 8.59480E-01]
Σₜᵢ_M87 = [2.16280E-01, 3.61700E-01, 5.05630E-01, 6.11170E-01, 5.08900E-01, 9.26670E-01, 9.60990E-01]
Σₜᵢ_UO2 = [2.12450E-01, 3.55470E-01, 4.85540E-01, 5.59400E-01, 3.18030E-01, 4.01460E-01, 5.70610E-01]
Σₜᵢ_Mod = [2.30070E-01, 7.76460E-01, 1.48420E+00, 1.50520E+00, 1.55920E+00, 2.02540E+00, 3.30570E+00]

frac_of_MFP = 1//4
materials[1].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_FC )) # Fission Chamber
materials[2].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_GT )) # Guide Tube
materials[3].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_M43)) # MOX-4.3% enriched
materials[4].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_M70)) # MOX-7.0% enriched
materials[5].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_M87)) # MOX-8.7% enriched
materials[6].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_UO2)) # Uranium Oxide
materials[7].mesh_size = frac_of_MFP*inv(mean(Σₜᵢ_Mod)) # Moderator

#gmsh.model.mesh.set_size(gmsh.model.get_entities(0), 0.05)
#gmsh.model.mesh.generate(2)
set_mesh_field_using_materials(materials)
#println("Mesh->Define->Size fields. Then, new view of Min field.")
#gmsh.fltk.run()
##gmsh.view.write(0, filename[1:end-5]*".pos")
#field_data = gmsh.view.get_model_data(0, 0)
generate_mesh(order = 1, faces = "Triangle", opt_iters = 2)
gmsh.write("c5g7.inp")
mesh_error = get_cad_to_mesh_error()
for i in eachindex(mesh_error)
    println(mesh_error[i])
end
#gmsh.fltk.run()
gmsh.finalize()

mesh = import_mesh("c5g7.inp")
statistics(mesh)
# Partition mesh according to mpact grid hierarchy, and write as an xdmf file
mpt = MeshPartitionTree(mesh)
export_mesh(mpt, "c5g7.xdmf")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
const ACTIVE_POWER = "P"
const ACTIVE_POWER_IN = "Pin"
const ACTIVE_POWER_OUT = "Pout"
const AREA_MISMATCH = "area_mismatch"
const COLD_START = "start_cold"
const ENERGY = "E"
const ENERGY_BUDGET = "energy_budget"
const FLOW_ACTIVE_POWER = "Fp"
const HOT_START = "start_hot"
const INFLOW = "In"
const ON = "On"
const REACTIVE_POWER = "Q"
const RESERVE = "R"
const SERVICE_REQUIREMENT = "service_requirement"
const SLACK_DN = "γ⁻"
const SLACK_UP = "γ⁺"
const SPILLAGE = "Sp"
const START = "start"
const STOP = "stop"
const THETA = "theta"
const VM = "Vm"
const WARM_START = "start_warm"
const LIFT = "z"

abstract type VariableType end

"""
Struct to dispatch the creation of Active Power Variables
"""
struct ActivePowerVariable <: VariableType end

function make_variable_name(
    ::Type{ActivePowerVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "P")
end

"""
Struct to dispatch the creation of Active Power Input Variables for 2-directional devices.
for instance storage or pump-hydro
"""
struct ActivePowerInVariable <: VariableType end

function make_variable_name(
    ::Type{ActivePowerInVariable},
    ::Type{T},
) where {T <: PSY.Device}
    return encode_symbol(T, "Pin")
end

"""
Struct to dispatch the creation of Active Power Output Variables for 2-directional devices.
for instance storage or pump-hydro
"""
struct ActivePowerOutVariable <: VariableType end

function make_variable_name(
    ::Type{ActivePowerOutVariable},
    ::Type{T},
) where {T <: PSY.Device}
    return encode_symbol(T, "Pout")
end

struct ColdStartVariable <: VariableType end

function make_variable_name(::Type{ColdStartVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "start_cold")
end

struct EnergyVariable <: VariableType end

function make_variable_name(::Type{EnergyVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "E")
end

struct HotStartVariable <: VariableType end

function make_variable_name(::Type{HotStartVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "start_hot")
end

struct LiftVariable <: VariableType end

function make_variable_name(::Type{LiftVariable})
    return encode_symbol("lift")
end

function make_variable_name(::Type{LiftVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "lift")
end

struct OnVariable <: VariableType end

function make_variable_name(::Type{OnVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "On")
end

struct ReactivePowerVariable <: VariableType end

function make_variable_name(
    ::Type{ReactivePowerVariable},
    ::Type{T},
) where {T <: PSY.Device}
    return encode_symbol(T, "Q")
end

struct ReserveVariable <: VariableType end

function make_variable_name(::Type{ReserveVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "R")
end

struct ActiveServiceVariable <: VariableType end

struct ServiceRequirementVariable <: VariableType end

function make_variable_name(
    ::Type{ServiceRequirementVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "service_requirement")
end

struct SpillageVariable <: VariableType end

function make_variable_name(::Type{SpillageVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "Sp")
end

struct StartVariable <: VariableType end

function make_variable_name(::Type{StartVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "start")
end

struct StopVariable <: VariableType end

function make_variable_name(::Type{StopVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, "stop")
end

struct WarmStartVariable <: VariableType end

function make_variable_name(::Type{WarmStartVariable}, ::Type{T}) where {T <: PSY.Device}
    return encode_symbol(T, WARM_START)
end

##### AGC Variables #####
struct SteadyStateFrequencyDeviation <: VariableType end

function make_variable_name(::Type{SteadyStateFrequencyDeviation})
    return encode_symbol("Δf")
end

struct AreaMismatchVariable <: VariableType end

function make_variable_name(::Type{AreaMismatchVariable})
    return encode_symbol(AREA_MISMATCH)
end

struct DeltaActivePowerUpVariable <: VariableType end

function make_variable_name(
    ::Type{DeltaActivePowerUpVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "ΔP_up")
end

function make_variable_name(
    ::Type{DeltaActivePowerUpVariable},
    ::Type{PSY.RegulationDevice{T}},
) where {T <: PSY.Device}
    return encode_symbol(T, "ΔP_up")
end

struct DeltaActivePowerDownVariable <: VariableType end

function make_variable_name(
    ::Type{DeltaActivePowerDownVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "ΔP_dn")
end

function make_variable_name(
    ::Type{DeltaActivePowerDownVariable},
    ::Type{PSY.RegulationDevice{T}},
) where {T <: PSY.Device}
    return encode_symbol(T, "ΔP_dn")
end

struct AdditionalDeltaActivePowerUpVariable <: VariableType end

function make_variable_name(
    ::Type{AdditionalDeltaActivePowerUpVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "ΔPe_up")
end

function make_variable_name(
    ::Type{AdditionalDeltaActivePowerUpVariable},
    ::Type{PSY.RegulationDevice{T}},
) where {T <: PSY.Device}
    return encode_symbol(T, "ΔPe_up")
end

struct AdditionalDeltaActivePowerDownVariable <: VariableType end

function make_variable_name(
    ::Type{AdditionalDeltaActivePowerDownVariable},
    ::Type{T},
) where {T <: PSY.Component}
    return encode_symbol(T, "ΔPe_dn")
end

function make_variable_name(
    ::Type{AdditionalDeltaActivePowerDownVariable},
    ::Type{PSY.RegulationDevice{T}},
) where {T <: PSY.Device}
    return encode_symbol(T, "ΔPe_dn")
end

struct SmoothACE <: VariableType end

function make_variable_name(
    ::Type{SmoothACE},
    ::Type{T},
) where {T <: PSY.AggregationTopology}
    return encode_symbol(T, "SACE")
end

function encode_symbol(::Type{T}, name1::AbstractString, name2::AbstractString) where {T}
    return Symbol(join((name1, name2, IS.strip_module_name(T)), PSI_NAME_DELIMITER))
end

function encode_symbol(
    ::Type{T},
    name1::AbstractString,
    name2::AbstractString,
) where {T <: PSY.Reserve}
    T_ = replace(IS.strip_module_name(T), "{" => "_")
    T_ = replace(T_, "}" => "")
    return Symbol(join((name1, name2, T_), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name1::Symbol, name2::Symbol) where {T}
    return encode_symbol(IS.strip_module_name(T), string(name1), string(name2))
end

function encode_symbol(::Type{T}, name::AbstractString) where {T}
    return Symbol(join((name, IS.strip_module_name(T)), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name::AbstractString) where {T <: PSY.Reserve}
    T_ = replace(IS.strip_module_name(T), "{" => "_")
    T_ = replace(T_, "}" => "")
    return Symbol(join((name, T_), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name::Symbol) where {T}
    return encode_symbol(T, string(name))
end

function encode_symbol(name::AbstractString)
    return Symbol(name)
end

function encode_symbol(name1::AbstractString, name2::AbstractString)
    return Symbol(join((name1, name2), PSI_NAME_DELIMITER))
end

function encode_symbol(name::Symbol)
    return name
end

function decode_symbol(name::Symbol)
    return split(String(name), PSI_NAME_DELIMITER)
end

make_variable_name(var_type, device_type) = encode_symbol(device_type, var_type)
make_variable_name(var_type) = encode_symbol(var_type)


text:
using ProgressBars

"""
    launch_experiment!(;
        ValueSelectionArray::Array{ValueSelection, 1},
        problem_type::Symbol=:coloring,
        problem_params::Dict=coloring_params,
        nbEpisodes::Int64=10,
        strategy::Type{DFSearch}=DFSearch,
        variableHeuristic=selectVariable
)

This functions launch an amount of nbEpisodes problems solving. The problems are created by
the given generator. The strategy used during the CP Search and the variable heuristic used can
be precised as well. To finish with, the value selection heuristic (eather learned or basic) are
given to function. Each problem generated will be solved once by every value selection heuristic
given, making it possible to compare them.

All the results of consecutive search during the training is stocked in metricsArray containing a metrics
for each heuristic (eather learned or basic).

It is also possible to give an evaluator to compare the evolution of performance of the heuristic on same instances along the search.

This function is called by `train!` and by `benchmark_solving!`.

Every "evalFreq" episodes, all heuristic are evaluated ( weights are no longer updated during the evaluation).
"""
function launch_experiment!(
        valueSelectionArray::Array{T, 1},
        generator::AbstractModelGenerator,
        nbEpisodes::Int64,
        strategy::S1,
        eval_strategy::S2,
        variableHeuristic::AbstractVariableSelection,
        out_solver::Bool,
        verbose::Bool;
        metrics::Union{Nothing, AbstractMetrics}=nothing,
        evaluator::Union{Nothing, AbstractEvaluator}=SameInstancesEvaluator(valueSelectionArray,generator),
        restartPerInstances::Int64,
    ) where{T <: ValueSelection, S1,S2 <: SearchStrategy}

    nbHeuristics = length(valueSelectionArray)

     #get the type of CPmodel ( does it contains an objective )
    trailer = Trailer()
    model = CPModel(trailer)
    fill_with_generator!(model, generator)
    metricsArray=AbstractMetrics[]
    for j in 1:nbHeuristics
        if !isnothing(metrics)
            push!(metricsArray,metrics(model,valueSelectionArray[j]))
        else
            push!(metricsArray,BasicMetrics(model,valueSelectionArray[j]))
        end
    end

    iter = ProgressBar(1:nbEpisodes)
    for i in iter
    #for i in 1:nbEpisodes
        verbose && println(" --- EPISODE: ", i)

        empty!(model)

        fill_with_generator!(model, generator)

        for j in 1:nbHeuristics
            reset_model!(model)
            if isa(valueSelectionArray[j], LearnedHeuristic)
                verbose && print("Visited nodes with learnedHeuristic : " )
            else
                verbose && print("Visited nodes with basic Heuristic n°$(j-1) : ")
            end
            dt = @elapsed for k in 1:restartPerInstances
                restart_search!(model)
                search!(model, strategy, variableHeuristic, valueSelectionArray[j], out_solver=out_solver)

                verbose && print(model.statistics.numberOfNodesBeforeRestart, ", ")
            end
            metricsArray[j](model,dt)  #adding results in the metrics data structure
            println()
        end

        if !isnothing(evaluator) && (i % evaluator.evalFreq == 0)
            evaluate(evaluator, variableHeuristic, eval_strategy; verbose = verbose)
        end
        verbose && println()
    end

    if !isnothing(evaluator)
        return metricsArray, evaluator.metrics

    end
    return metricsArray,[]
end


text:
# import package
julia> using Pkg

# package status
julia> Pkg.status()
    Status `C:\Users\praveen\.julia\environments\v1.0\Project.toml`
  [7876af07] Example v0.5.1

# add CSV package
julia> Pkg.add("CSV")
  Updating registry at `C:\Users\praveen\.julia\registries\General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
    Fetching: [===>                                     ]  6.3 %

# alternately
# or you can use the Pkg REPL mode by pressing ] (similar to ?, help mode, and ;, shell mode):

(v1.0) pkg> add "CSV"

using CSV, DataFrames, BenchmarkTools
iris = CSV.read(joinpath(dirname(pathof(DataFrames)),"..","test/data/iris.csv"))


text:
# push!(LOAD_PATH, "/home/caseykneale/Desktop/ChemometricsTools/ChemometricsTools.jl/");
# using Pkg
# Pkg.activate(".")
using Documenter, ChemometricsTools

makedocs(
	sitename = "ChemometricsTools",
	modules = [ChemometricsTools],
	format = Documenter.HTML(# Use clean URLs, unless built as a "local" build
	        prettyurls = !("local" in ARGS),
		),
	authors = "Casey Kneale",
    doctest = true,
	pages = ["Home" => "index.md",
			"Demos" => Any[
						"Transforms" => "Demos/Transforms.md",
						"Pipelines" => "Demos/Pipelines.md",
						"Classification" => "Demos/ClassificationExample.md",
						"Regression" => "Demos/RegressionExample.md",
						"SIPLS" => "Demos/SIPLS.md",
						"Calibration Transfer" => "Demos/CalibXfer.md",
						"Curve Resolution" => "Demos/CurveResolution.md"
					   ],
			 "Data" => Any[
			 					"Data Utilities" => "man/datautils.md",
			 			 		"Kernel Density Generator" => "man/kerneldensity.md"
							],
			 "Manipulations" => Any[
			 					"Preprocessing" => "man/Preprocess.md",
			 			 		"Transformations/Pipelines" => "man/Transformations.md",
						 		"Sampling" => "man/Sampling.md"
							],
			 "Visualizations" => Any[ "Plotting" => "man/Plotting.md"],
			 "Analysis" => Any[
			 					"DOE" => "man/DOE.md",
								"Univariate" => "man/Univariate.md",
			 					"Distances/Kernels" => "man/Dists.md",
			 					"Analysis" => "man/Analysis.md",
			  					"Clustering" => "man/Clustering.md",
								"Stats." => "man/Stats.md"
							],
			 "Modeling" => Any[
			 					"Training" => "man/Training.md",
						 		"Time Series" => "man/TimeSeries.md",
						 		"Regression Models" => "man/RegressionModels.md",
						 		"Regression Metrics" => "man/regressMetrics.md",
			 			 		"Classification Models" => "man/ClassificationModels.md",
						 		"Classification Metrics" => "man/classMetrics.md",
						 		"Tree Methods" => "man/Trees.md",
						 		"Ensemble Models" => "man/Ensemble.md"
							],
			 "Speciality Tools" => Any[
			 					 "Model Analysis" => "man/modelanaly.md",
						 		 "MultiWay" => "man/MultiWay.md",
								 "Hyperspectral" => "man/Hyperspectral.md",
								 "Anomaly Detection" => "man/AnomalyDetection.md",
								 "Curve Resolution" => "man/CurveResolution.md"
							],
			"Advanced" => Any[
								 "PSO" => "man/PSO.md",
								 "Genetic Algorithms" => "man/GeneticAlgorithms.md"
		    			]
		]
)

deploydocs(
    repo = "github.com/caseykneale/ChemometricsTools.jl.git",
	branch = "gh-pages",
)


text:
module NeatText

# Export for Removers
export remove_patterns,remove_stopwords,fix_contractions,remove_terms_in_bracket
export remove_puncts,remove_numbers,remove_phonenumbers,remove_currencies,remove_currency_symbols,remove_special_characters
export remove_emails,remove_urls,remove_hashtags,remove_htmltags,remove_userhandles,remove_emojis
export remove_md5sha,remove_btc,remove_pobox,remove_mastercards,remove_visacards,remove_streetaddress

# Export for Removers!(Mutating Fxns)
export remove_patterns!,remove_stopwords!,fix_contractions!,remove_terms_in_bracket!
export remove_puncts!,remove_numbers!,remove_phonenumbers!,remove_currencies!,remove_currency_symbols!,remove_special_characters!
export remove_emails!,remove_urls!,remove_hashtags!,remove_htmltags!,remove_userhandles!,remove_emojis!
export remove_md5sha!,remove_btc!,remove_pobox!,remove_mastercards!,remove_visacards!,remove_streetaddress!

# Export for Extractors
export extract_patterns,extract_stopwords
export extract_puncts,extract_numbers,extract_phonenumbers,extract_currencies,extract_dates
export extract_emails,extract_urls,extract_hashtags,extract_htmltags,extract_emojis,extract_userhandles
export extract_mastercards,extract_visacards,extract_creditcards,extract_md5sha,extract_pobox,extract_streetaddress
export extract_terms_in_bracket

# Export for Misc Functions
export clean_text!,clean_text

# Constants
const PUNCTS_REGEX = Regex("[!&'()*,-./:;?@[\\]^_`{|}]")
const EMAIL_REGEX =  r"[\w\.-]+@[\w\.-]+"
const NUMBERS_REGEX = r"\d+"
const PHONENUMBERS_REGEX = r"[\+\(]?[1-9][0-9 .\-\(\)]{8,}[0-9]"
const SPECIAL_CHARACTERS_REGEX = r"[^A-Za-z0-9 ]+"
const DATE_REGEX = r"([0-9]{2}\/[0-9]{2}\/[0-9]{4})|([0-9]{4}\/[0-9]{2}\/[0-9]{2})"

# modified for full stops, question marks, commas, colons and semi-colons, exclamation marks and quotation marks.
const MOST_COMMON_PUNCT_REGEX = Regex("""[!"&',-.;?_`]""")
const HTML_TAGS_REGEX = r"<[^<]+?>"
const CURRENCY_REGEX = r"[$¢£¤¥ƒ֏؋৲৳૱௹฿៛ℳ元円圆圓﷼\u20A0-\u20C0]\d+"
const CURRENCY_SYMB_REGEX = r"[$¢£¤¥ƒ֏؋৲৳૱௹฿៛ℳ元円圆圓﷼\u20A0-\u20C0]"
const HASHTAG_REGEX = r"#\S+"
const USER_HANDLES_REGEX = r"@\S+"
const STREET_ADDRESS_REGEX  = r"\d{1,4} [\w\s]{1,20}(?:street|st|avenue|ave|road|rd|highway|hwy|square|sq|trail|trl|drive|dr|court|ct|park|parkway|pkwy|circle|cir|boulevard|blvd)\W?`"
const PoBOX_REGEX  = r"(?i)P\.? ?O\.? Box \d+"
const ZIP_REGEX  = r"\b\d{5}(?:[-\s]\d{4})?\b"
const BTC_ADDRESS_REGEX = r"[13][a-km-zA-HJ-NP-Z1-9]{25,34}"
const CreditCard_REGEX = r"(?:(?:(?:\d{4}[- ]?){3}\d{4}|\d{15,16}))"
const VISACard_REGEX = r"4\d{3}[\s-]?\d{4}[\s-]?\d{4}[\s-]?\d{4}"
const MASTERCard_REGEX  = r"5[1-5]\d{2}[\s-]?\d{4}[\s-]?\d{4}[\s-]?\d{4}"
const MD5_SHA_REGEX = r"([0-9a-fA-F]{32})|([0-9a-fA-F]{40})|([0-9a-fA-F]{64})"
# const URL_REGEX= r"(?:^|(?<![\w\/\.]))|(?:(?:https?:\/\/|ftp:\/\/|www\d{0,3}\.))"
const URL_REGEX = r"http[s]?://(?:[a-zA-Z]|[0-9]|[$-_@.&+]|[!*\(\),]|(?:%[0-9a-fA-F][0-9a-fA-F]))+"

const EMOJI_REGEX = Regex("""[
                       u"\U0001F600-\U0001F64F"  # for emoticons
                       u"\U0001F300-\U0001F5FF"  # for symbols & pictographs
                       u"\U0001F680-\U0001F6FF"  # for transport & map symbols
                       u"\U0001F1E0-\U0001F1FF"  # for flags (iOS)
                       u"\U00002702-\U000027B0"
                       u"\U000024C2-\U0001F251"
                       "]+""",)

"""
remove_patterns: remove from given text your specified patterns
params:
	-text: text or string of text to be cleaned

"""
function remove_patterns(text::AbstractString,pattern::Regex)
	return replace(text,pattern=>" ")
end

"""
remove_puncts: Remove punctuations from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_puncts(text::AbstractString)
	return replace.(text,PUNCTS_REGEX=>" ")
end

"""
remove_special_characters: Remove Special Characters from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_special_characters(text::AbstractString)
	return replace.(text,SPECIAL_CHARACTERS_REGEX=>" ")
end

"""
remove_emails: Remove emails from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_emails(text::AbstractString)
	return replace.(text,EMAIL_REGEX=>" ")
end

"""
remove_emojis: Remove emojis from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_emojis(text::AbstractString)
	return replace.(text,EMOJI_REGEX=>" ")
end

"""
remove_numbers: Remove numbers from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_numbers(text::AbstractString)
	return replace.(text,NUMBERS_REGEX=>" ")
end

"""
remove_phonenumbers: Remove phonenumbers from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_phonenumbers(text::AbstractString)
	return replace.(text,PHONENUMBERS_REGEX=>" ")
end

"""
remove_htmltags: Remove htmltags from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_htmltags(text::AbstractString)
	return replace.(text,HTML_TAGS_REGEX=>" ")
end

"""
remove_urls: Remove URLS from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_urls(text::AbstractString)
	return replace.(text,URL_REGEX=>" ")
end

"""
remove_dates: Remove Dates from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_dates(text::AbstractString)
	return replace.(text,DATE_REGEX=>" ")
end

"""
remove_currencies: Remove currencies from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_currencies(text::AbstractString)
	return replace.(text,CURRENCY_REGEX=>" ")
end

"""
remove_currency_symbols: Remove currency symbols from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_currency_symbols(text::AbstractString)
	return replace.(text,CURRENCY_SYMB_REGEX=>" ")
end

"""
remove_userhandles: Remove userhandles or mentions from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_userhandles(text::AbstractString)
	return replace.(text,USER_HANDLES_REGEX=>" ")
end

"""
remove_hashtags: Remove hashtags symbols from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_hashtags(text::AbstractString)
	return replace.(text,HASHTAG_REGEX=>" ")
end

"""
remove_streetaddress: Remove street addresses from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_streetaddress(text::AbstractString)
	return replace.(text,STREET_ADDRESS_REGEX=>" ")
end

"""
remove_pobox: Remove Post Office Box Address from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_pobox(text::AbstractString)
	return replace.(text,PoBOX_REGEX=>" ")
end

"""
remove_btc: Remove BIT Address from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_btc(text::AbstractString)
	return replace.(text,BTC_ADDRESS_REGEX=>" ")
end

"""
remove_mastercard: Remove Master Card from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_mastercards(text::AbstractString)
	return replace.(text,MASTERCard_REGEX=>" ")
end

"""
remove_visacard: Remove Visa Card numbers from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_visacards(text::AbstractString)
	return replace.(text,VISACard_REGEX=>" ")
end

"""
remove_creditcards: Remove CreditCard numbers from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_creditcards(text::AbstractString)
	return replace.(text,CreditCard_REGEX=>" ")
end

"""
remove_md5sha: Remove MD5 and SHA values from a given text
params:
	-text: text or string of text to be cleaned

"""
function remove_md5sha(text::AbstractString)
	return replace.(text,MD5_SHA_REGEX =>" ")
end

"""
fix_contractions: fix contractions for a given text

"""
function fix_contractions(text::AbstractString)
	text = replace.(text,"ain't"=>"am not")
	text = replace.(text,"aren't"=>"are not")
	text = replace.(text,"can't"=>"cannot")
	text = replace.(text,"can't've"=>"cannot have")
	text = replace.(text,"'cause"=>"because")
	text = replace.(text,"could've"=>"could have")
	text = replace.(text,"couldn't"=>"could not")
	text = replace.(text,"couldn't've"=>"could not have")
	text = replace.(text,"didn't"=>"did not")
	text = replace.(text,"doesn't"=>"does not")
	text = replace.(text,"don't"=>"do not")
	text = replace.(text,"hadn't"=>"had not")
	text = replace.(text,"hadn't've"=>"had not have")
	text = replace.(text,"hasn't"=>"has not")
	text = replace.(text,"haven't"=>"have not")
	text = replace.(text,"he'd"=>"he would")
	text = replace.(text,"he'd've"=>"he would have")
	text = replace.(text,"he'll"=>"he will")
	text = replace.(text,"he'll've"=>"he will have")
	text = replace.(text,"he's"=>"he is")
	text = replace.(text,"how'd"=>"how did")
	text = replace.(text,"how'd'y"=>"how do you")
	text = replace.(text,"how'll"=>"how will")
	text = replace.(text,"how's"=>"how is")
	text = replace.(text,"I'd"=>"I would")
	text = replace.(text,"I'd've"=>"I would have")
	text = replace.(text,"I'll"=>"I will")
	text = replace.(text,"I'll've"=>"I will have")
	text = replace.(text,"I'm"=>"I am")
	text = replace.(text,"I've"=>"I have")
	text = replace.(text,"isn't"=>"is not")
	text = replace.(text,"it'd"=>"it had")
	text = replace.(text,"it'd've"=>"it would have")
	text = replace.(text,"it'll"=>"it will")
	text = replace.(text,"it'll've"=>"it will have")
	text = replace.(text,"it's"=>"it is")
	text = replace.(text,"let's"=>"let us")
	text = replace.(text,"ma'am"=>"madam")
	text = replace.(text,"mayn't"=>"may not")
	text = replace.(text,"might've"=>"might have")
	text = replace.(text,"mightn't"=>"might not")
	text = replace.(text,"mightn't've"=>"might not have")
	text = replace.(text,"must've"=>"must have")
	text = replace.(text,"mustn't"=>"must not")
	text = replace.(text,"mustn't've"=>"must not have")
	text = replace.(text,"needn't"=>"need not")
	text = replace.(text,"needn't've"=>"need not have")
	text = replace.(text,"o'clock"=>"of the clock")
	text = replace.(text,"oughtn't"=>"ought not")
	text = replace.(text,"oughtn't've"=>"ought not have")
	text = replace.(text,"shan't"=>"shall not")
	text = replace.(text,"sha'n't"=>"shall not")
	text = replace.(text,"shan't've"=>"shall not have")
	text = replace.(text,"she'd"=>"she would")
	text = replace.(text,"she'd've"=>"she would have")
	text = replace.(text,"she'll"=>"she will")
	text = replace.(text,"she'll've"=>"she will have")
	text = replace.(text,"she's"=>"she is")
	text = replace.(text,"should've"=>"should have")
	text = replace.(text,"shouldn't"=>"should not")
	text = replace.(text,"shouldn't've"=>"should not have")
	text = replace.(text,"so've"=>"so have")
	text = replace.(text,"so's"=>"so is")
	text = replace.(text,"that'd"=>"that would")
	text = replace.(text,"that'd've"=>"that would have")
	text = replace.(text,"that's"=>"that is")
	text = replace.(text,"there'd"=>"there had")
	text = replace.(text,"there'd've"=>"there would have")
	text = replace.(text,"there's"=>"there is")
	text = replace.(text,"they'd"=>"they would")
	text = replace.(text,"they'd've"=>"they would have")
	text = replace.(text,"they'll"=>"they will")
	text = replace.(text,"they'll've"=>"they will have")
	text = replace.(text,"they're"=>"they are")
	text = replace.(text,"they've"=>"they have")
	text = replace.(text,"to've"=>"to have")
	text = replace.(text,"wasn't"=>"was not")
	text = replace.(text,"we'd"=>"we had")
	text = replace.(text,"we'd've"=>"we would have")
	text = replace.(text,"we'll"=>"we will")
	text = replace.(text,"we'll've"=>"we will have")
	text = replace.(text,"we're"=>"we are")
	text = replace.(text,"we've"=>"we have")
	text = replace.(text,"weren't"=>"were not")
	text = replace.(text,"what'll"=>"what will")
	text = replace.(text,"what'll've"=>"what will have")
	text = replace.(text,"what're"=>"what are")
	text = replace.(text,"what's"=>"what is")
	text = replace.(text,"what've"=>"what have")
	text = replace.(text,"when's"=>"when is")
	text = replace.(text,"when've"=>"when have")
	text = replace.(text,"where'd"=>"where did")
	text = replace.(text,"where's"=>"where is")
	text = replace.(text,"where've"=>"where have")
	text = replace.(text,"who'll"=>"who will")
	text = replace.(text,"who'll've"=>"who will have")
	text = replace.(text,"who's"=>"who is")
	text = replace.(text,"who've"=>"who have")
	text = replace.(text,"why's"=>"why is")
	text = replace.(text,"why've"=>"why have")
	text = replace.(text,"will've"=>"will have")
	text = replace.(text,"won't"=>"will not")
	text = replace.(text,"won't've"=>"will not have")
	text = replace.(text,"would've"=>"would have")
	text = replace.(text,"wouldn't"=>"would not")
	text = replace.(text,"wouldn't've"=>"would not have")
	text = replace.(text,"y'all"=>"you all")
	text = replace.(text,"y'alls"=>"you alls")
	text = replace.(text,"y'all'd"=>"you all would")
	text = replace.(text,"y'all'd've"=>"you all would have")
	text = replace.(text,"y'all're"=>"you all are")
	text = replace.(text,"y'all've"=>"you all have")
	text = replace.(text,"you'd"=>"you had")
	text = replace.(text,"you'd've"=>"you would have")
	text = replace.(text,"you'll"=>"you you will")
	text = replace.(text,"you'll've"=>"you you will have")
	text = replace.(text,"you're"=>"you are")
	text = replace.(text,"you've"=>"you have")
	return text
end

"""Returns a string with terms/words in bracket together with the bracket removed
Params:
    - bracket_form: string or symbol about whether is it  curly({) or square([])
            or etc
"""

function remove_terms_in_bracket(text::AbstractString;bracket_form="curly")
        CURLY_BRACKET_REGEX = r"\{(.*?)\}"
        SQUARE_BRACKET_REGEX = r"\[(.*?)\]"
        NORMAL_BRACKET_REGEX = r"\((.*?)\)"
        if (bracket_form == "curly") | (bracket_form == "{}")
            result = replace.(text,CURLY_BRACKET_REGEX=>"")
        elseif (bracket_form == "square") | (bracket_form == "[]")
            result = replace.(text,SQUARE_BRACKET_REGEX=>"")
        elseif (bracket_form == "normal") | (bracket_form == "()")
            result = replace.(text,NORMAL_BRACKET_REGEX=>"")
           end
        return result
end

function count_char(s::AbstractString)
  res = Dict{Char, Int}()
 for c in s
      if haskey(res, c)
         res[c] += 1
      else
         res[c] = 1
      end
   end
   return res
end

# A Dictionary of Regular Expressions
const REGEX_DICT = Dict(
	:puncts => PUNCTS_REGEX,
	:special_char => SPECIAL_CHARACTERS_REGEX,
	:userhandles => USER_HANDLES_REGEX,
	:emails => EMAIL_REGEX,
	:numbers => NUMBERS_REGEX,
	:phonenumbers => PHONENUMBERS_REGEX,
	:mentions => r"@\w+",
	:hashtags => r"#\w+",
	:urls => r"http[s]?://(?:[a-zA-Z]|[0-9]|[$-_@.&+]|[!*\(\),]|(?:%[0-9a-fA-F][0-9a-fA-F]))+",
	:htmltags => HTML_TAGS_REGEX
);

"""
clean_text: Clean A Given Text using a set of patterns

usepatterns: An array or collection of predefined pattern symbols
	- [:emails,:urls,:mentions:hashtags,:userhandles,:htmltags:puncts,:numbers,:phonenumbers,:special_char]

   [:emails
	:urls
	:mentions
	:hashtags
	:userhandles
	:htmltags
	:puncts
	:numbers
	:phonenumbers
	:special_char]

note:
----
	- the order of the patterns matters in how the text is cleaned
"""
function clean_text(text::AbstractString;usepatterns::Array)
	for p in usepatterns
		text = replace.(text,REGEX_DICT[p]=> "")
	end
	return text
end

include("stopwordlist.jl")

"""
remove_stopwords: returns a text with all stopwords in the given text removed
param:
-----
 - lang: specify the language for the stopwords : en,es,fr,ru,yo,de
params:
	-text: text or string of text to be cleaned

"""
function remove_stopwords(text::AbstractString,lang="en")
	if lang == "en"
		stopwords_in_use = STOPWORDS_en
    elseif lang == "es"
    	stopwords_in_use = STOPWORDS_es
    elseif lang == "fr"
    	stopwords_in_use = STOPWORDS_fr
    elseif lang == "ru"
    	stopwords_in_use = STOPWORDS_ru
    elseif lang == "yo"
    	stopwords_in_use = STOPWORDS_yo
    elseif lang == "de"
    	stopwords_in_use = STOPWORDS_de
    else
    	stopwords_in_use = STOPWORDS_en
    end

	alltokens = []
	for w in split(lowercase(text)," ")
		if w ∉ stopwords_in_use
			push!(alltokens,w)
		end
	end
	return join(alltokens," ")

end

# Extractors
include("textExtractors.jl")

end # module


text:
using MakieGraphs
using Test

using GLMakie
using GraphPlot
using LightGraphs

println("nothing here yet")


text:
# # Root Mean Squared Fluctuation (RMSF)
#
# md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__cookbook/notebooks/03_RMSF.ipynb)
# md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__cookbook/notebooks/03_RMSF.ipynb)
#
#
# ## Problem description
#
# The [Root Mean Squared Fluctuation (RMSF)](https://en.wikipedia.org/wiki/Mean_squared_displacement)
# is a common way to measure residue flexibility in a structural ensemble.
# It is a measure of how far is the residue moving from its average position
# in the group of structures. Usually, we represent a residue position with
# the spatial coordinates of its alpha carbon.
#
# The protein structures should be previously superimposed to calculate the
# RMSF, for example, by using the `superimpose` function of the
# [`PDB` module of `MIToS`](@ref Module-PDB). In this example, we are going
# to measure the RMSF of each residue from an NMR ensemble using the
# `rmsf` function.
#
# The structure superimposition could be the most complicated step of the
# process, depending on the input data. In particular, it structures come
# from different PDB structures or homologous proteins can require the use
# of external programs,
# as [MAMMOTH-mult](https://ub.cbm.uam.es/software/online/mamothmult.php) or
# [MUSTANG](https://lcb.infotech.monash.edu/mustang/) among others,
# tailored for this task.
#
# In this case, we are going to use an NMR ensemble. Therefore, we are not
# going to need to superimpose the structures as NMR models have the
# same protein sequence and are, usually, well-aligned.
#
#
# ## MIToS solution

import MIToS
using MIToS.PDB
using Plots

# Lets read the NMR ensemble:

pdb_file   = abspath(pathof(MIToS), "..", "..", "test", "data", "1AS5.pdb")
pdb_res = read(pdb_file, PDBFile, occupancyfilter=true)
# md nothing # hide

# We set `occupancyfilter` to `true` to ensure that we have one single set of
# coordinates for each atom. That filter isn't essential for NMR structures,
# but It can avoid multiple alpha carbons in crystallographic structures with
# disordered atoms.
# We can get an idea of the alpha carbon positions by plotting these residues:

scatter(pdb_res, legend=false)

# As we saw in the previous plot, the structure doesn't need to be
# superimposed. Now, we are going to separate each model into different
# vectors, storing each vector into a `Dict`:

models = Dict{String,Vector{PDBResidue}}()
for res in pdb_res
	push!(get!(models, res.id.model, []), res)
end

# Then, we simply need to collect all the PDB models in the values
# of the `Dict`, to get the vector of `PDBResidues` vectors required
# to calculate the RMSF.

pdb_models = collect(values(models))
# md nothing # hide

# And, finally, call the `rmsf` function on the list of structures. It is
# important that all the vectors has the same number of `PDBResidue`s. This
# function assumes that the nth element of each vector corresponds to the same
# residue:

RMSF = rmsf(pdb_models)

# This return the vector of RMSF values for each residue, calculated using
# the coordinates of the alpha carbons.
# You can plot this vector to get an idea of the which are the most flexible
# position in your structure:

plot(RMSF, legend=false, xlab="Residue", ylab="RMSF [Å]")


text:
using ArgParse2
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using BioAlignments
using DataFramesMeta
using Statistics
using st16SseqJuliaTools
using ProgressMeter
using XAM

function main()
    parser = ArgumentParser(prog = "analyse_alnself.jl",
                        description = "Detect chimeras",
                        epilog = "Removing potential chimeras and repair snps")
    add_argument!(parser, "--ref", "-r", help = "reference fasta",type = String)
    add_argument!(parser, "--ideal", "-i", help = "ideal contigs - for debugging",type = String,default="")
    add_argument!(parser, "--salmon-rez", "-s", help = "Salmon rezults",type = String)
    add_argument!(parser, "--out", "-o", help = "Fasta file to write out cleaned up sequences",default="chimerasfree.fasta",type = String)
    #add_argument!(parser, "--output", "-o", help = "File with names",type = String)
    args = parse_args(parser)
    salmon_data = read_salmon(args.salmon_rez)
    ref = fasta_to_dict(args.ref)
    length_data=Dict{String,Int64}()
    for k in keys(ref)
        length_data[k] = length(ref[k])
    end

    #Limits for real match ... mayve should be options
    minl = 100
    minid = 0.96
    pretrim_length = 250
    minid_for_close = 0.96 # Identity of first 250 fragments of contigs
    max_unaligned_length_to_consider_conservative_start_extention = 20
    min_dif_in_expression = 0.1
    min_dif_in_length = 0.1

    costmodel = BioAlignments.CostModel(match=0, mismatch=1, insertion=1, deletion=1);
    #known_chimeras = ["1_1868;size=92", "1_6093;size=54", "1_2188;size=22","1_1760;size=8"]
    args.ideal == "" ? known_good_contigs = [] : known_good_contigs = collect(keys(fasta_to_dict(args.ideal)))
    pairs_4chim = Array{Tuple{String,String,Int64,Int64,Float64,Bool},1}()
    pairs_4contain = Array{Tuple{String,String,Int64,Int64,Float64,Float64,Int64},1}()
    aln_costs = AffineGapScoreModel(
        match=5,
        mismatch=-4,
        gap_open=-10,
        gap_extend=-5
    )

    ref_keys = collect(keys(ref))
    @showprogress 1 "Pairwise alignments..."  for k1 in ref_keys
        #TO DO - below there are three consequtive alignments with some similar steps - should be a function there
        pairs_part_4chim = Array{Tuple{String,String,Int64,Int64,Float64,Bool},1}(undef,length(ref_keys))
        pairs_part_4contain = Array{Tuple{String,String,Int64,Int64,Float64,Float64,Int64},1}(undef,length(ref_keys))
        Threads.@threads for i in 1:length(ref_keys)
            k2 = ref_keys[i]
            if k1 != k2
                rn = k1
                tn = tempname = k2
                rs = ref[rn]
                ts = ref[tn]
                ts_wopref = ts[pretrim_length+1:end]
                aln = alignment(pairalign(LocalAlignment(), rs, ts_wopref, aln_costs))
                matches = BioAlignments.count_matches(aln)
                mismatches = BioAlignments.count_mismatches(aln)
                ins = count_insertions(aln)
                del = count_deletions(aln)
                s = pos = first(aln.a.aln.anchors).seqpos+1
                e = rpos = last(aln.a.aln.anchors).seqpos
                l = abs(s-e)+1
                idfrac = matches/l
                t_pos=first(aln.a.aln.anchors).refpos+1
                t_rpos=last(aln.a.aln.anchors).refpos
                tseq = ref[tempname][pretrim_length+1:end][t_pos:t_rpos]
                rtseq = rs[pos:rpos]
                #check if  close sequences and the detected match could be due to extention of conservative region
                t_pref = ref[tempname][1:250]
                r_pref = ref[rn][1:250]
                pairaln_pref = BioAlignments.pairalign(BioAlignments.EditDistance(),t_pref,r_pref,costmodel)
                pairaln_pref = BioAlignments.alignment(pairaln_pref)
                matches_pref = BioAlignments.count_matches(pairaln_pref)
                idfrac_pref = matches_pref/250
                is_cons5region_extention = (idfrac_pref >= minid_for_close && abs(t_pos+pretrim_length-pos)<= max_unaligned_length_to_consider_conservative_start_extention)
                # Calculate overall alignment
                aln_fulllen = alignment(pairalign(LocalAlignment(), rs, ts, aln_costs))
                matches_fulllen = BioAlignments.count_matches(aln_fulllen)
                mismatches_fulllen = BioAlignments.count_mismatches(aln_fulllen)
                ins_fulllen = count_insertions(aln_fulllen)
                del_fulllen = count_deletions(aln_fulllen)
                s_fulllen = pos_fulllen = first(aln_fulllen.a.aln.anchors).seqpos
                e_fulllen = rpos_fulllen = last(aln_fulllen.a.aln.anchors).seqpos
                t_pos_fulllen=first(aln_fulllen.a.aln.anchors).refpos+1
                t_rpos_fulllen=last(aln_fulllen.a.aln.anchors).refpos
                taln_fulllen = abs(t_pos_fulllen-t_rpos_fulllen)+1
                idfrac_fulllen = matches_fulllen/taln_fulllen

                if  (rn in ["5_20;size=1"]) && (tempname in["5_3133;size=1"])
                     println()
                     println(aln)
                     println()
                     println(aln_fulllen)
                     println(pos_fulllen," ",rpos_fulllen," ",length(rs)," ",length(ts))
                     println((tempname,rn,t_pos_fulllen,t_rpos_fulllen,idfrac_fulllen,taln_fulllen/length(ts),length(rs)))
                end
                if taln_fulllen > 0
                    pairs_part_4contain[i] = (tempname,rn,t_pos_fulllen,t_rpos_fulllen,idfrac_fulllen,taln_fulllen/length(ts),length(rs))
                end
                if idfrac >= minid &&
                    l >= minl
                    if  (rn in ["1_2540;size=1"]) && (tempname in["6_421;size=1"])
                        println("UUUUUUU ",is_cons5region_extention)
                        println(idfrac_pref,"   ",idfrac," AAAAAAA")
                        println("********************************************")
                        println("PAIR: \n>r\n$rtseq\n>t\n$tseq\n---------------")
                        println(">tf")
                        println(ref[tempname][pretrim_length+1:end])
                        println(">rf")
                        println(ref[rn])
                        println(">tf_250")
                        println(ref[tempname][1:250])
                        println(">rf_250")
                        println(ref[rn][1:250])
                        println(">tf_full_$tempname")
                        println(ref[tempname])
                        println(">rf_full_$rn")
                        println(ref[rn])
                        println(tempname," ",rn," ",t_pos," ",t_rpos," ",pos," ",rpos)
                        println("AAAAAAAAAAAAAAAAAAAaaa",idfrac)
                        println(aln)
                        println("aln of prefix")
                        println(pairaln_pref)
                        println(idfrac_pref)
                        println(minid_for_close)
                        println(t_pos+pretrim_length-pos)
                        println(max_unaligned_length_to_consider_conservative_start_extention)
                    end
                    pairs_part_4chim[i] = (tempname,rn,t_pos,t_rpos,idfrac,is_cons5region_extention)
                end
            end
        end
        for ii in 1:length(ref_keys)
            if isassigned(pairs_part_4chim,ii)
                push!(pairs_4chim,pairs_part_4chim[ii])
            end
            if isassigned(pairs_part_4contain,ii)
                push!(pairs_4contain,pairs_part_4contain[ii])
            end
        end
    end
    #Collect data for  chimera detection
    tn_rn_map = Dict{String,Array{String,1}}()
    for pair in pairs_4chim
        tn,rn,is_cons5region_extention =  pair[1],pair[2],pair[6]
        if !is_cons5region_extention
            if !(tn in keys(tn_rn_map))
                tn_rn_map[tn] = Array{String,1}()
            end
            push!(tn_rn_map[tn],rn)
        end
    end
    #detect if chimeric or not
    chimeras = Array{String,1}()
    #println("--------------------------------------------------------------")
    #testf = open("test.csv","w")
    #println(testf,"tn max_r_id  good_t good_r good_both chimeric false_chimeric  tn_s max_r_s difs len_t len_r difl  ")
    for tn in keys(tn_rn_map)
        tn_s = salmon_data[split(tn,";")[1]]
        rn_ss = Array{Float64,1}()
        for rn in tn_rn_map[tn]
            rn_s = salmon_data[split(rn,";")[1]]
            push!(rn_ss, rn_s)
        end
        pot_chim = DataFrame(ref_names = tn_rn_map[tn],ref_salmon = rn_ss)
        sort!(pot_chim,[:ref_salmon],rev=true)
        max_r_s = pot_chim[1,:].ref_salmon
        max_r_id = pot_chim[1,:].ref_names
        len_t = length_data[tn]
        len_r = length_data[max_r_id]
        difl = (len_t-len_r)/len_t

        chimeric = (tn_s < max_r_s) && (len_t <= len_r)

        good_t = (tn in known_good_contigs )
        good_r = (max_r_id in known_good_contigs )
        good_both = good_t && good_r
        difs = abs(max_r_s - tn_s )/tn_s
        false_chimeric = chimeric && good_t
        ct = 0

       # println(testf,"$tn $max_r_id  $good_t $good_r $good_both $chimeric $false_chimeric  $tn_s $max_r_s $difs $len_t $len_r $difl  ")
        if  chimeric
            println("Template: ID: $tn  TPM: $tn_s  good: $good_t  Most abundant chimeric partner:  ID: $max_r_id   TPM: $max_r_s   good: $good_r")
            ct += 1
            push!(chimeras,tn)
            if chimeric && good_t
                println("##############MISCLASSIFICATION##############")
                println("
                Template:
                    ID: $tn
                    TPM: $tn_s
                    good: $good_t
                Most abundant chimeric partner:
                    ID: $max_r_id
                    TPM: $max_r_s
                    good: $good_r
                ")
                ts = ref[tn]
                rs = ref[max_r_id]
                pairalnres = BioAlignments.pairalign(BioAlignments.EditDistance(),rs,ts[pretrim_length+1:end],aln_costs)
                pairaln = BioAlignments.alignment(pairalnres)
                println(pairaln)
                println("Exiting...")
                exit()
            end
        end

    end
    #close(testf)

    chimeras = unique(chimeras)

    #Collect data for  contained detection

    contained = Array{String,1}()
    tn_rn_map = Dict{String,Array{Tuple{String,String,Float64,Float64,Int64},1}}()
    for pair in pairs_4contain
        tn,rn,idf,covf,rl =  pair[1],pair[2],pair[5],pair[6],pair[7]
        if !(tn in keys(tn_rn_map))
            tn_rn_map[tn] = Array{Tuple{String,String,Float64,Float64,Int64},1}()
        end
        push!(tn_rn_map[tn],(tn,rn,idf,covf,rl))
    end

    #detect contained
    #merge perfect matches per high coverage
    contained_idf_limit_c = 0.96 # Limid for identity for containment
    contained_covf_limit_c = 0.99 # Limid for coverage for containment
    #merge perfect matches per high identity fragment - merge ideal matches allowing sme artifacts at the ends...
    contained_idf_limit_i = 0.99 # Limid for identity for containment
    contained_covf_limit_i = 0.95 # Limid for coverage for containment
    for ts in tn_rn_map["5_3133;size=1"]
        println(ts)
    end
    for tn in keys(tn_rn_map)
        tl = length(ref[tn])
        cont = false
        for rd in tn_rn_map[tn]
            if tl <= rd[5]
                if tn == "5_3133;size=1i" && rd[2] == "5_20;size=1"
                    println("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA")
                    println((round(rd[3],digits=2) >= contained_idf_limit_c && round(rd[4],digits=2) >= contained_covf_limit_c))
                    println((round(rd[3],digits=2), contained_idf_limit_i, round(rd[4],digits=2) ,contained_covf_limit_i))
                    exit()
                end
                if (round(rd[3],digits=2) >= contained_idf_limit_c && round(rd[4],digits=2) >= contained_covf_limit_c) ||
                (round(rd[3],digits=2) >= contained_idf_limit_i && round(rd[4],digits=2) >= contained_covf_limit_i)
                    if tl == rd[5]
                        #if lengths matches add abundance criteria
                        r_s = salmon_data[split(rd[2],";")[1]]
                        t_s = salmon_data[split(tn,";")[1]]
                        if r_s > t_s
                            cont = true
                            break
                        end
                    else
                        cont = true
                        break
                    end
                end
            end
        end
        if cont
            push!(contained,tn)
        end
    end
    filtered_out = OrderedDict{String,LongSequence{BioSequences.DNAAlphabet{4}}}()
    ct_chimeric = 0
    ct_contained = 0
    ct_good = 0
    for k in keys(ref)
        ischim = (k in chimeras)
        iscont = (k in contained)
        if k == "5_3133;size=1"
            println("AAAAAAAAAAAH ",ischim, " ",iscont)
        end
        if !ischim && !iscont
            ct_good += 1
            filtered_out[k]= ref[k]
        else
            if ischim
                ct_chimeric += 1
            end
            if iscont
                ct_contained += 1
            end
        end
    end
    dict_to_fasta(filtered_out,args.out)
    println("$ct_chimeric chimeras , $ct_contained contained sequences were detected, $ct_good sequences were  written to ",args.out)
end

main()


text:
using LightGBM

@testset "stringifyparams -- convert to zero-based" begin
    indices = [1, 3, 5, 7, 9]
    classifier = LightGBM.LGBMClassification(categorical_feature = indices)
    ds_parameters = LightGBM.stringifyparams(classifier, LightGBM.DATASETPARAMS)

    expected = "categorical_feature=0,2,4,6,8"
    @test occursin(expected, ds_parameters)
end


text:
### LLL method

### f = prod(fs) over Zp, find factors over Z
function mu(i,j,B,Bs)
    dot(B[i,:], Bs[j,:]) // dot(Bs[j,:], Bs[j,:] )
end
function swaprows!(A, i, j, tmp)
    tmp[:] = A[i,:]
    A[i,:] = A[j,:]
    A[j,:] = tmp
    nothing
end

function _reduce(i, T, B, U)
    j = i - 1
    while j > 0
        rij = round(T, digits=U[i,j])
        B[i,:] -= rij * B[j, :]
        U[i,:] -= rij * U[j, :]
        j -= 1
    end
end

# Find a short vector for B = [f1;f2;...;fn]
# short vector if f1 after the work is done
# http://algo.epfl.ch/_media/en/projects/bachelor_semester/rapportetiennehelfer.pdf .
# Algo 2
function LLL_basis_reduction!(B::Matrix{T}, c=4//3) where {T <: Integer} # c > 4/3

    m, n = size(B)

    Bstar =  Matrix{Rational{T}}(I, m, m)
    U = Matrix{Rational{T}}(I, m, m)

    tmpZ = zeros(T, 1, m)
    tmpQ = zeros(Rational{T}, 1, m)

    ONE, ZERO = one(Rational{T}), zero(Rational{T})

    ## initialize Bstar, U
    Bstar[1,:] = B[1,:]  # b^*_1 = b_1
    for i in 2:n
        Bstar[i,:] = B[i,:]
        for j in 1:(i-1)
            U[i,j] = mu(i,j,B,Bstar) #vecdot(bi, bstar_j) // vecdot(bstar_j, bstar_j)
            Bstar[i,:] -= U[i,j] * Bstar[j,:]
        end

        _reduce(i, T, B, U)
    end

    i = 1
    while i < m
        if norm(B[i,:],2) < c * norm(B[i+1,:],2)
            i += 1
        else
            # Modify Q and R in order to keep the relation B = QR after the swapping
            Bstar[i+1, :] += U[i+1,i] * Bstar[i, :]

            U[i, i] = mu(i, i+1, B, Bstar)
            U[i, i+1] = U[i+1, i] = ONE
            U[i+1, i+1] = ZERO

            Bstar[i, :] -= U[i,i] * Bstar[i+1, :]

            swaprows!(U, i, i+1, tmpQ)
            swaprows!(Bstar, i, i+1, tmpQ)
            swaprows!(B, i, i+1, tmpZ)

            for k in (i+2):m
                U[k,i] = mu(k, i, B, Bstar)
                U[k,i+1] = mu(k, i+1, B, Bstar)
            end

            if abs(U[i+1,i]) > 1//2
                _reduce(i+1, T, B, U)
            end
            i = max(i-1, 1)
        end
    end
end

## For a polynomial p, create a lattice of p, xp, x^2p, ... then find a short vector.
## This will be a potential factor of f
function short_vector(u::AbstractAlgebra.Generic.Poly, m, j, d)
    us = poly_coeffs(u)
    T = eltype(us)
    F = zeros(T, j, j)

    # coefficients of us*x^i
    for i in 0:(j-d-1)
        F[1 + i, (1 + i):(1+i+d)] = us
    end
    # coefficients of mx^i
    for i in 0:(d-1)
        F[1 + j - d + i, 1 + i] = m
    end

    LLL_basis_reduction!(F)
    x = variable(u)
    # short vector is first row
    as_poly(vec(F[1,:]), x)
end

## u is a factor over Zp
## we identify irreducible gstar of f
## return gstar, fstar = f / gstar (basically), and reduced Ts
function identify_factor(u::AbstractAlgebra.Generic.Poly, fstar::AbstractAlgebra.Generic.Poly, Ts::Vector, p::T, l, b, B) where {T}

    d, nstar = degree(u), degree(fstar)
    d < nstar || throw(DomainError())
    m = p^l
    j = d + 1
    gstar = u
    Tsp = as_poly_Zp.(Ts, p)

    while j <= nstar+1
        gstar = short_vector(u,m,j,d)

        gstarp = as_poly_Zp(gstar, p)
        inds = Int[]
        for (i,t) in enumerate(Tsp)
            divides(gstarp, Tsp[i])[1] && push!(inds, i)  # allocates
        end

        notinds = setdiff(1:length(Ts), inds)
        Ssp = Tsp[notinds]

        # find hstar = prod(Ts)
        hstar = as_poly(mod(b, m) * (isempty(Ssp) ? one(gstarp) : prod(Ssp)))

        pd = onenorm(primpart(gstar)) * onenorm(primpart(hstar))

        if pd < B
            gstar, hstar = primpart(gstar), primpart(hstar)
            Ts = Ts[notinds] # modify Ts (Ts[:])or make a copy?
            fstar = primpart(hstar)
            break
        end
        j = j + 1
    end

    primpart(gstar), fstar, Ts
end

"""
Use LLL basis reduction to identify factors of square free f from its factors over Zp
"""
function identify_factors_lll(f, facs, p, l, b, B)

    Ts = sort(facs, by=degree, rev=true)

    Gs = similar(Ts, 0)
    fstar = primpart(f)
    m = p^l

    while length(Ts) > 0
        u = popfirst!(Ts)
        if length(Ts) == 0
            push!(Gs, fstar)
            break
        end
        if degree(fstar) == degree(u)
             push!(Gs, fstar)
            break
        end
        #        return (u, fstar, Ts, p, l, b, B)
        gstar, fstar, Ts = identify_factor(u, fstar, Ts, p, l, b, B)
        push!(Gs, gstar)
    end
    Gs
end


text:
function rfloor(val::Float64)::Integer
    rf_val = Integer(floor(val + val * 1e-10 + 1e-6))
    if (rf_val < val - 1 + 1e-6)
        rf_val += 1
    end
    return rf_val
end

function scale_vector_of_double(v::Vector{Float64})
    n = length(v)
    int_v = Vector{Integer}()
    max_val = maximum(v)
    scaling_factor = typemax(Cint) / (n + 1) / max_val
    for i in 1:n
        push!(int_v, rfloor(scaling_factor * v[i]))
    end
    return int_v
end

function minknap(p::Vector{Float64}, w::Vector{T}, capacity::T) where T <: Integer
    int_p = scale_vector_of_double(p)
    (obj, sol) = minknap(int_p, w, capacity)
    real_obj = sum(sol[i] * p[i] for i in 1:length(sol))
    return (real_obj, sol)
end

function minmcknap(p::Vector{Float64}, w::Vector{T}, capacity::T, items_per_class::Vector{T}) where T <: Integer
    int_p = scale_vector_of_double(p)
    (obj, sol) = minmcknap(int_p, w, capacity, items_per_class)
    real_obj = sum(sol[i] * p[i] for i in 1:length(sol))
    return (real_obj, sol)
end