text:
import JuMP

@testset "shapetrait and nin/nout" begin
    import NaiveNASlib: shapetrait, NoShapeTraitError, noutnin_errmsg
    struct NoutNinFailDummy end

    @test shapetrait(NoutNinFailDummy()) == NoutNinFailDummy

    @test_throws NoShapeTraitError nout(NoutNinFailDummy())
    @test_throws NoShapeTraitError nin(NoutNinFailDummy())

    @test noutnin_errmsg(NoutNinFailDummy(), nout) == "nout not defined for NoutNinFailDummy()! Either implement NaiveNASlib.nout(::NoutNinFailDummy) or implement NaiveNASlib.nout(::ST, ::NoutNinFailDummy) where ST is the return type of NaiveNASlib.shapetrait(::NoutNinFailDummy)"
    @test noutnin_errmsg(NoutNinFailDummy(), nin) == "nin not defined for NoutNinFailDummy()! Either implement NaiveNASlib.nin(::NoutNinFailDummy) or implement NaiveNASlib.nin(::ST, ::NoutNinFailDummy) where ST is the return type of NaiveNASlib.shapetrait(::NoutNinFailDummy)"

    struct NinNoutSuccessDummy end
    NaiveNASlib.nout(::NinNoutSuccessDummy, args...) = 3
    NaiveNASlib.nin(::NinNoutSuccessDummy, args...) = [3]

    @test nout(NinNoutSuccessDummy()) == 3
    @test nin(NinNoutSuccessDummy()) == [3]

    # This is not really how the traits are supposed to be used, but it saves us from defining a new type and associated nout and nin methods here
    struct NinNoutTraitSuccessDummy end
    NaiveNASlib.shapetrait(::NinNoutTraitSuccessDummy) = NinNoutSuccessDummy()

    @test nout(NinNoutTraitSuccessDummy()) == 3
    @test nin(NinNoutTraitSuccessDummy()) == [3]
end

@testset "Size mutations" begin

    inpt(size, id="in") = inputvertex(id, size)
    nt(name) = named(name)
    tf(name) = t -> nt(name)(AfterΔSizeTrait(validateafterΔsize(base=ΔSizeFailNoOp()), t))
    cc(in...; name="cc") = conc(in..., dims=1, traitdecoration = tf(name))
    ea(ins...; name="ea") = +(traitconf(tf(name)) >> ins[1], ins[2:end]...)

    @testset "Basic tests" begin
        using NaiveNASlib: Input, Output, Both, neighbours

        av(in, outsize, name="av") = absorbvertex(SizeDummy(nout(in), outsize), in; traitdecoration=tf(name))

        @testset "AbsorbVertex" begin
            v0 = av(inpt(2), 2, "iv")
            v1 = av(v0, 3, "v1")

            @test nin(v1) == [2]
            @test nout(v1) == 3

            @test outputs.(inputs(v1)) == [[v1]]

            Δnin!(v1, 2)
            @test nin(v1) == [4]

            Δnin!(v1, -3)
            @test nin(v1) == [1]

            Δnout!(v1, -2)
            @test nout(v1) == 1

            # Add one vertex and see that change propagates
            v2 = av(v1, 4, "v2")
            @test nout(v1) == nin(v2)[1]
            @test outputs(v1) == [v2]
            @test inputs(v2) == [v1]

            Δnin!(v2, 4)
            @test nout(v1) == nin(v2)[1] == 5

            Δnout!(v1, -2)
            @test nout(v1) == nin(v2)[1] == 3

            # Fork of v1 into a new vertex
            v3 = av(v1, 2, "v3")
            @test outputs(v1) == [v2, v3]
            @test inputs(v3) == inputs(v2) == [v1]

            Δnout!(v1, -2)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 1

            Δnin!(v3, 3)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 4

            Δnin!(v2, -2)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 2
        end

        @testset "StackingVertex" begin

            @testset "StackingVertex 1 to 1" begin

                v0 = av(inpt(2), 2, "v0")
                v1 = cc(v0)
                v2 = av(v1, 2,"v2")

                @test outputs.(inputs(v2)) == [[v2]]
                @test outputs.(inputs(v1)) == [[v1]]
                @test outputs(v0) == [v1]
                @test neighbours(Both(), v1) == [v0, v2]

                Δnout!(v0, 3)
                @test [nout(v0)] == nin(v1) == [nout(v1)] == nin(v2) == [5]

                Δnin!(v2, -2)
                @test [nout(v0)] == nin(v1) == [nout(v1)] == nin(v2) == [3]

                Δnin!(v1, +1)
                @test [nout(v0)] == nin(v1) == [nout(v1)] == nin(v2) == [4]

                Δnout!(v1, -1)
                @test [nout(v0)] == nin(v1) == [nout(v1)] == nin(v2) == [3]
            end

            @testset "StackingVertex 2 inputs" begin
                # Try with two inputs to StackingVertex
                v1 = av(inpt(2,"in1"), 2, "v1")
                v2 = av(inpt(3,"in2"), 3, "v2")
                v3 = cc(v1, v2)
                v4 = av(v3, 5, "v4")

                @test inputs(v3) == [v1, v2]
                @test outputs.(inputs(v3)) == [[v3], [v3]]
                @test outputs(v1) == [v3]
                @test outputs(v2) == [v3]
                @test neighbours(Both(), v3) == [v1, v2, v4]

                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == [5]

                Δnout!(v2, -2)
                @test nin(v4) == [nout(v3)] == [sum(nin(v3))] == [nout(v1) + nout(v2)]  ==  [3]

                Δnin!(v4, 3)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == [6]

                Δnout!(v1, 4)
                Δnin!(v4, -8)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == [2]
                @test nout(v1) == nout(v2) == 1

                #Add another output
                v5 = av(v3, 2, "v5")

                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [2]

                Δnout!(v1, 3)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [5]

                Δnin!(v5, -2)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [3]

                Δnin!(v4, 3)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [6]

                Δnin!(v3, -1, missing)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [5]

                Δnin!(v3, missing, 2)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [7]

                Δnout!(v3, -1)
                @test [nout(v1) + nout(v2)] == [sum(nin(v3))] == [nout(v3)] == nin(v4) == nin(v5) == [6]
            end

        end
        @testset "InvariantVertex" begin

            @testset "InvariantVertex 1 to 1" begin
                #Behaviour is identical to StackingVertex in this case
                v1 = av(inpt(2), 2, "v1")
                v2 = ea(v1)
                v3 = av(v2, 2, "v3")

                @test outputs.(inputs(v3)) == [[v3]]
                @test outputs.(inputs(v2)) == [[v2]]
                @test outputs(v1) == [v2]

                Δnout!(v1, 3)
                @test [nout(v1)] == nin(v2) == [nout(v2)] == nin(v3) == [5]

                Δnin!(v3, -2)
                @test [nout(v1)] == nin(v2) == [nout(v2)] == nin(v3) == [3]
            end

            @testset "invariantvertex 2 inputs" begin
                # Try with two inputs a v3ariantVertex
                v1 = av(inpt(2, "in1"), 3, "v1")
                v2 = av(inpt(2, "in2"), 3, "v2")
                v3 = ea(v1, v2)
                v4 = av(v3, 5, "v4")

                @test inputs(v3) == [v1, v2]
                @test outputs.(inputs(v3)) == [[v3], [v3]]
                @test outputs(v1) == [v3]
                @test outputs(v2) == [v3]

                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == 3

                Δnout!(v2, -2)
                @test nin(v4)[] == nout(v3) == nin(v3)[1] == nin(v3)[2]  == nout(v1) == nout(v2)  ==  1

                Δnin!(v4, 3)
                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == 4

                #Add another output
                v5 = av(v3, 4, "v5")

                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == nin(v5)[] == 4

                Δnout!(v1, -3)
                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == nin(v5)[] == 1

                Δnin!(v5, 2)
                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == nin(v5)[] == 3

                Δnout!(v2, 3)
                @test nout(v1) == nout(v2) == nin(v3)[1] == nin(v3)[2] == nout(v3) == nin(v4)[] == nin(v5)[] == 6
            end

            @testset "invariantvertex with params and size" begin
                struct InvariantWithSize
                    s::Int
                end
                NaiveNASlib.nin(i::InvariantWithSize, ::SizeInvariant, v::AbstractVertex) = [i.s]

                v = invariantvertex(InvariantWithSize(7), inpt(3, "in"); traitdecoration=tf("v"))

                @test nin(v) == [nout(v)] == [7]
            end
        end

        @testset "validate after Δsize" begin
            import NaiveNASlib: ΔSizeFailError

            # Mock for creating invalid size changes
            struct Ignore{T}
                f::T
            end
            NaiveNASlib.nout(i::Ignore) = nout(i.f)
            NaiveNASlib.nin(i::Ignore) = nin(i.f)
            NaiveNASlib.Δsize!(::Ignore, ::AbstractVector{<:Integer}, ::Integer) = nothing
            NaiveNASlib.Δsizetype(i::Ignore) = NaiveNASlib.Δsizetype(i.f)

            @testset "With AfterΔSizeTrait" begin

                v1 = inpt(3, "v1")
                v2 = absorbvertex(Ignore(SizeDummy(nout(v1), 5)), v1; traitdecoration=tf("v2"))
                v3 = av(v2, 7, "v3")
                v4 = absorbvertex(Ignore(SizeDummy(nout(v3), 5)), v3; traitdecoration=tf("v4"))

                @test_throws ΔSizeFailError Δnout!(v2 => 2)
                @test_throws ΔSizeFailError Δnin!(v3 => -3)
                @test_throws ΔSizeFailError Δnin!(v4 => 5)
                @test_throws ΔSizeFailError Δnout!(v3 => -7)
            end

            @testset "As strategy" begin
                v1 = inpt(3, "v1")
                v2 = absorbvertex(Ignore(SizeDummy(nout(v1), 5)), v1; traitdecoration=nt("v2"))
                v3 = absorbvertex(SizeDummy(nout(v2), 7), v2; traitdecoration=nt("v3"))
                v4 = absorbvertex(Ignore(SizeDummy(nout(v3), 5)), v3; traitdecoration=nt("v4"))

                @test_throws ΔSizeFailError Δsize!(validateafterΔsize(base=ΔNout(v2 => 2)))
                @test_throws ΔSizeFailError Δsize!(validateafterΔsize(base=ΔNin(v3 =>-3)))
                @test_throws ΔSizeFailError Δsize!(validateafterΔsize(base=ΔNin(v4, 5)))
                @test_throws ΔSizeFailError Δsize!(validateafterΔsize(base=ΔNout(v3 => -7)))
            end
        end

        @testset "Fail mutate immutable" begin
            import NaiveNASlib: ΔSizeFailError

            v0 = inpt(3)
            v1 = av(v0, 5, "v1")

            @test_throws ΔSizeFailError Δsize!(ΔNoutExact(v0, 2;fallback=ThrowΔSizeFailError("Success?!")), [v0, v1])
            @test [nout(v0)] == nin(v1) == [3]

            @test_throws ΔSizeFailError Δsize!(ΔNinExact(v1, -2; fallback=ThrowΔSizeFailError("Success?!")), [v1, v0])
            @test [nout(v0)] == nin(v1) == [3]
        end
    end

    @testset "Time limit" begin
        import NaiveNASlib: TimeLimitΔSizeStrategy, sizemodel
        @test JuMP.time_limit_sec(sizemodel(TimeLimitΔSizeStrategy(123), AbstractVertex[])) == 123
    end

    @testset "Mutate tricky structures" begin

        ## Helper functions
        iv(in; name="iv") = invariantvertex(identity, in; traitdecoration=tf(name))
        concd1(paths...;name="conc") = conc(paths..., dims=1, traitdecoration=tf(name))
        av(in, outsize, name="comp") = absorbvertex(SizeDummy(nout(in), outsize), in, traitdecoration = tf(name))
        function stack(start, nouts...; bname = "stack")
            # Can be done on one line with mapfoldl, but it is not pretty...
            next = start
            for i in 1:length(nouts)
                next = av(next, nouts[i], "$(bname)_$i")
            end
            return next
        end

        @testset "Residual fork block" begin
            start = av(inputvertex("in", 3), 9, "start")
            p1 = stack(start, 3,4, bname = "p1")
            p2 = stack(start, 4,5, bname = "p2")
            resout = "add" >> start + concd1(p1, p2)
            out = av(resout, 9, "out")

            @test nout(resout) == 9

            # Propagates to out, start outputs of start and also to p1 as objective is minimized by increasing p1 by 1
            @test Δnout!(p2, -2)
            @test nout(p2) == 3
            @test nin(out) == [nout(start)] == [8]
            @test nout(p2) + nout(p1) == 8
            #outputs(start) = first vertex in p1, p2 and resout (which has two inputs)
            @test reduce(vcat, nin.(outputs(start))) == [8, 8, 8, 8]

            @test Δnin!(out, 2)
            @test nin(out) == [nout(start)] == [10]
            @test reduce(vcat, nin.(outputs(start))) == [10, 10, 10, 10]
        end

        @testset "Half transparent residual fork block" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 4, "split")
            p1 = iv(split, name="p1") #Just an identity vertex
            p2 = stack(split, 3,2,4, bname="p2")
            resout = "add" >> start + concd1(p1, p2)
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Propagates to input of first vertex of p2, input of out and start
            # via p1 and resout as well as to input of split
            @test Δnout!(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [8]
            @test reduce(vcat, nin.(outputs(split))) == [3, 3]

            # Should basically undo the previous mutation
            @test Δnin!(out, +1)
            @test nin(out) == [nout(start)] == nin(split) == [9]
        end

        @testset "Transparent fork block" begin
            start = av(inputvertex("in", 3), 4, "start")
            p1 = iv(start, name="p1")
            p2 = iv(start, name="p2")
            joined = concd1(p1, p2, name="join")
            out = av(joined, 3, "out")

            @test nout(joined) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout.
            @test Δnout!(start, -1)
            @test nin(out) == [2nout(start)] == [6]

            # Should basically undo the previous mutation
            @test Δnin!(out, 2)
            @test nin(out) == [2nout(start)] == [8]
        end

        @testset "Transparent residual fork block" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 4, "split")
            p1 = iv(split, name="p1")
            p2 = iv(split, name="p2")
            resout = "add" >> start + concd1(p1, p2, name="join")
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout.
            @test Δnout!(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [6]

            # Should basically undo the previous mutation
            @test Δnin!(out, +2)
            @test nin(out) == [nout(start)] == nin(split) == [8]

            @test Δnout!(start, -2)
            @test nin(out) == [nout(start)] == [6]
            @test nout(split) == 3
        end

        @testset "Transparent residual fork block with single absorbing path" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 3, "split")
            p1 = iv(split, name="p1")
            p2 = iv(split, name="p2")
            p3 = av(split, 2, "p3")
            resout = "add" >> start + concd1(p1, p2, p3, name="join")
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout.
            @test Δnout!(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [6]

            # Should basically undo the previous mutation
            @test Δnin!(out, 2)
            @test nin(out) == [nout(start)] == nin(split) == [8]

            @test Δnout!(start, -2)
            @test nin(out) == [nout(start)] == [2*nout(split) + nout(p3)] == [6]
            @test nout(split) == 2
            @test nout(p3) == 2
        end
    end

    @testset "Size Mutation possibilities" begin
        using NaiveNASlib: newsizes

        # Helpers
        struct SizeConstraint{T}
            constraint::Int
            w::T
        end

        function NaiveNASlib.compconstraint!(::NaiveNASlib.ScalarSize, s, c::SizeConstraint, data)
            fv_out = JuMP.@variable(data.model, integer=true)
            JuMP.@constraint(data.model, c.constraint * fv_out ==  nout(data.vertex) - data.noutdict[data.vertex])

            ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))
            fv_in = JuMP.@variable(data.model, [1:length(ins)], integer=true)
            JuMP.@constraint(data.model, [i=1:length(ins)], c.constraint * fv_in[i] ==  nout(ins[i]) - data.noutdict[ins[i]])
        end
        NaiveNASlib.Δsizetype(c::SizeConstraint) = NaiveNASlib.Δsizetype(c.w)
        NaiveNASlib.nout(c::SizeConstraint) = nout(c.w)
        NaiveNASlib.nin(c::SizeConstraint) = nin(c.w)
        NaiveNASlib.Δsize!(c::SizeConstraint, insize::AbstractVector{<:Integer}, outsize::Integer) = Δsize!(c.w, insize, outsize)
        av(size, csize, in... ;name = "av") = absorbvertex(SizeConstraint(csize, SizeDummy(collect(nout.(in)), size)), in..., traitdecoration=tf(name))

        @testset "SizeStack multi inputs" begin
            v1 = av(6,3, inpt(3), name="v1")
            v2 = av(7,2, inpt(3), name="v2")
            v3 = av(8,2, inpt(3), name="v3")

            cc1 = cc(v1, v2, name="cc1")
            cc2 = cc(v3, v2, v1, v2, name="cc2")

            # Expect only v1 to change as size change is not compatible with v2
            @test Δnout!(cc1, -3)
            @test nin(cc1) == [nout(v1), nout(v2)] == [3, 7]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [8, 7, 3, 7]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

            # v1 can't change as it is too small already
            @test Δnout!(cc2, -6)
            @test nin(cc1) == [nout(v1), nout(v2)] == [3, 5]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [6, 5, 3, 5]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

            # Evil action! Must have understanding that change will propagate to cc1 from
            # both v1 and v2
            @test Δnout!(cc2, 9)
            @test nin(cc1) == [nout(v1), nout(v2)] == [6, 7]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [8, 7, 6, 7]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

        end

        @testset "Stacked SizeStacks" begin
            v1 = av(100,1, inpt(3), name="v1")
            v2 = av(100,2, inpt(3), name="v2")
            v3 = av(100,3, inpt(3), name="v3")
            v4 = av(100,5, inpt(3), name="v4")

            cc1 = cc(v2, v3, name="cc1")
            cc2 = cc(cc1, v1, name="cc2")
            cc3 = cc(cc2, v4, v2, name="cc3")

            @test Δnout!(cc3, -60)
            @test nin(cc1) == nout.(inputs(cc1)) == [86, 91]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 87]
            @test nin(cc3) == nout.(inputs(cc3)) == [264, 90, 86]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 264
            @test nout(cc3) == nout(cc2) + nout(v4) + nout(v2) == 440

            v5 = av(10, 3, cc3, name="v5")

            @test Δnout!(v1, 3)
            @test nin(cc1) == nout.(inputs(cc1)) == [86, 91]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 90]
            @test nin(cc3) == nout.(inputs(cc3)) == [267, 90, 86]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 267
            @test [nout(cc3)] == [(nout(cc2) + nout(v4) + nout(v2))] == nin(v5) == [443]

            # Evil action! Must have understanding that the change in v2 will propagate
            # to cc3 input 1 through cc1 and cc2
            @test Δnout!(v2, -12)
            @test nin(cc1) == nout.(inputs(cc1)) == [74, 103]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 90]
            @test nin(cc3) == nout.(inputs(cc3)) == [267, 90, 74]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 267
            @test [nout(cc3)] == [(nout(cc2) + nout(v4) + nout(v2))] == nin(v5) == [431]
        end

        @testset "SizeInvariant multi input" begin
            v1 = av(10,1, inpt(3), name="v1")

            ea1 = ea(v1, name="ea1")
            ea2 = ea(ea1, v1, name="ea2")

            @test Δnout!(ea2, -2)

            @test nout(ea2) == nout(ea1) == nout(v1) == 8
            @test nin(ea2) == [nout(ea1), nout(v1)] == [8, 8]
            @test nin(ea1) == [nout(ea1)] == [8]
        end

        @testset "SizeInvariant multi SizeStack input" begin
            v1 = av(7,1, inpt(3), name="v1")
            v2 = av(13,1, inpt(3), name="v2")

            cc1 = cc(v1, v2, name="cc1")
            cc2 = cc(v1, v2, name="cc2")

            ea1 = ea(cc1, cc2, name="ea1")
            ea2 = ea(cc1, cc2, name="ea2")
            ea3 = ea(ea1, ea2, name="ea3")

            @test Δnout!(ea3, -4)

            @test nout(ea3) == nout(ea2) == nout(ea1) == nout(cc2) == nout(cc1) == nout(v2) + nout(v1) == 16
            @test nin(ea3) == nin(ea2) == nin(ea3) == [nout(cc1), nout(cc2)] == [16, 16]
            @test nin(cc2) == nin(cc1) == [nout(v1), nout(v2)] == [6, 10]
        end

        @testset "Stacked SizeInvariant" begin
            v1 = av(100,1, inpt(3), name="v1")
            v2 = av(100,2, inpt(3), name="v2")
            v3 = av(100,3, inpt(3), name="v3")
            v4 = av(100,5, inpt(3), name="v4")

            ea1 = ea(v2, v3, name="ea1")
            ea2 = ea(ea1, v1, name="ea2")
            ea3 = ea(ea2, v4, v2, name="ea3")

            # Everything touches everything in this setup
            @test Δnout!(ea3, -30)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 70
            @test nout(ea1) == nout(ea2) == nout(ea3) == 70

            v5 = av(10, 3, ea3, name="v5")

            @test Δnout!(v1, 60)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 130
            @test [nout(ea1)] == [nout(ea2)] == [nout(ea3)] == nin(v5) == [130]

            # Evil action! Must have understanding that the change in v2 will propagate
            # to ea3 input 1 through ea1 and ea2 and hold off updating it through input 3
            @test Δnout!(v2, -30)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 100
            @test [nout(ea1)] == [nout(ea2)] == [nout(ea3)] == nin(v5) == [100]
        end

        @testset "Stacked input to SizeInvariant" begin
            v0 = inpt(3, "in")
            v1 = av(20, 2, v0, name="v1")
            v2 = av(30, 3, v0, name="v2")
            v3 = cc(v1,v2,v1, name="v3")
            v4 = av(70, 5, v0, name="v4")
            v5 = ea(v4,v3, name="v5")
            v6 = av(10, 7, v5, name="v6")

            @test Δnout!(v5, 2*2*3*5*7)

            @test nout(v5) == nout(v4) == nout(v3) == 490
            @test nin(v5) == [nout(v4), nout(v3)] == [490, 490]
            @test nin(v3) == [nout(v1), nout(v2), nout(v1)] == [140, 210, 140]
        end

        @testset "Invariant input to SizeInvariant" begin
            v0 = inpt(3, "in")
            v1 = av(20, 2, v0, name="v1")
            v2 = av(20, 3, v0, name="v2")
            v3 = ea(v1,v2, name="v3")
            v4 = av(20, 5, v0, name="v4")
            v5 = ea(v4,v3, name="v5")
            v6 = av(10, 7, v5, name="v6")

            @test Δnout!(v5, 2*3*5*7)

            @test nout(v5) == nout(v4) == nout(v3) == 230
            @test nin(v5) == [nout(v4), nout(v3)] == nin(v3) == [nout(v1), nout(v2)] == [230, 230]
        end

        @testset "SizeStack duplicate through SizeInvariant" begin
            v0 = inpt(5, "in")
            v1 = av(2, 2, v0, name="v1")
            v2 = ea(v1, name="v2")
            v3 = ea(v1, name="v3")
            v4 = cc(v2,v3, name="v4")
            v5 = ea(v4, name="v5")
            v6 = cc(v5, v2, name="v6")

            @test Δnout!(v6, 12)

            @test nout(v6) == sum(nin(v6)) == 18
            @test nout(v5) == nin(v5)[] == 12
            @test nout(v4) == sum(nin(v4)) == 12
            @test nout(v1) == 6
        end

        @testset "SizeInvariant zig-zag" begin
                v0 = inpt(5, "in")
                v1 = av(2, 2, v0, name="v1")
                v2 = av(3, 3, v0, name="v2")
                v3 = cc(v1,v2, name="v3")
                function zigzag(vin1, sc, vin2=v0;name="zig")
                    vnew = av(nout(vin1), sc, vin2, name=name*"_new")
                    vout = ea(vnew, name=name*"_ivB")
                    vcon = ea(vout, vin1, name=name*"_ivA")
                    return vout
                end
                v4 = zigzag(v3, 5, name="z1")
                v5 = zigzag(v4, 7, name="z2")
                v6 = av(11, 11, v5, name="v6")

                expectedΔf = 2*3*5*7*11

                @test Δnout!(v2, expectedΔf)
                @test nout(v3) == nout(v4) == nout(v5) == expectedΔf + 2+3
                @test nin(v4) == nout.(inputs(v4)) == [expectedΔf + 2+3]
                @test nin(v5) == nout.(inputs(v5)) == [expectedΔf + 2+3]
        end

        @testset "SizeStack duplicate SizeInvariant mini-zig-zag" begin
            v0 = inpt(5, "in")
            v1 = av(2, 2, v0, name="v1")
            v2 = av(2, 3, v0, name="v2")
            v3 = ea(v1, name="v3")
            v4 = ea(v3, v2, name="v4")
            v5 = av(3, 5, v3, name="v5")
            v6 = cc(v4,v4, name="v6")
            v7 = ea(v6, v6, name="v7")
            v8 = av(4, 7, v7, name="v8")

            expectedΔf = 2*2*3*5*7

            @test Δnout!(v6, expectedΔf)
            @test nin(v8) == [nout(v7)] == [nout(v6)] == [4+expectedΔf]

            @test nin(v6) == [2 + expectedΔf ÷ 2, 2 + expectedΔf ÷ 2]
            @test nout(v1) == nout(v2) == 2 + expectedΔf ÷ 2
            @test nin(v5) == [2 + expectedΔf ÷ 2]
        end

        @testset "Deep SizeStack" begin
            v0 = inpt(3, "in")
            v1 = av(8, 1, v0, name="v1")
            v2 = av(4, 1, v0, name="v2")
            v3 = cc(v1,v2, name= "v3")
            pa1 = ea(v3, name="pa1")
            pb1 = ea(v3, name="pb1")
            pc1 = ea(v3, name="pc1")
            pd1 = av(5, 1, v3, name="pd1")
            pa1pa1 = ea(pa1, name="pa1pa1")
            pa1pb1 = ea(pa1, name="pa1pb1")
            pa2 = cc(pa1pa1, pa1pb1, name = "pa2")
            v4 = cc(pa2, pb1, pc1, pd1, name = "v4")

            @test Δnout!(v4, 8)

            @test nout(v4) == 61
            @test nin(v4) == nout.(inputs(v4)) == [28, 14, 14, 5]
            @test nin(pa2) == nout.(inputs(pa2)) == [14, 14]
            @test nin(pa1pa1) == nin(pa1pb1) == [14]
            @test nin(pa1) == nin(pb1) == nin(pc1) == [nout(v3)] == [14]
            @test nin(v3) == [10, 4]
        end

        @testset "Fail invalid size change" begin
            import NaiveNASlib: ΔSizeFailError
            v0 = inpt(3)
            v1 = av(3,3, v0, name="v1")
            v2 = av(5,2, v1, name="v2")

            @test_throws ΔSizeFailError newsizes(ΔNoutExact(v1, 2; fallback=ThrowΔSizeFailError("")), all_in_graph(v1))
            @test_throws ΔSizeFailError newsizes(ΔNoutExact(v1, 3; fallback=ThrowΔSizeFailError("")), all_in_graph(v1))

            @test_throws ΔSizeFailError newsizes(ΔNinExact(v2, 2; fallback=ThrowΔSizeFailError("")), all_in_graph(v1))
            @test_throws ΔSizeFailError newsizes(ΔNinExact(v2, 3; fallback=ThrowΔSizeFailError("")), all_in_graph(v1))

            @test newsizes(ΔNoutExact(v1, 2; fallback=ΔSizeFailNoOp()), all_in_graph(v1))[1] == false
            @test newsizes(ΔNinExact(v2, 2; fallback=ΔSizeFailNoOp()), all_in_graph(v2))[1] == false

            @test (@test_logs (:info, "Giving up") newsizes(ΔNoutExact(v1, 2; fallback=LogΔSizeExec("Giving up")), all_in_graph(v1)))[1] == false

            @test_logs (:warn, r"Could not change nout of .* by 2! Relaxing constraints...") Δnout!(v1, 2)
            @test [nout(v1)] == nin(v2) == [9]

            @test_logs (:warn, r"Could not change nin of .* by -2! Relaxing constraints...") Δnin!(v2, -2)
            @test [nout(v1)] == nin(v2) == [3]

            v3 = cc(v1,v2, name="v3")

            @test_logs (:warn, r"Could not change nin of .* by 2 and 6! Relaxing constraints...") Δnin!(v3, 2, 6)
            @test [nout(v1)] == nin(v2) == [9]
            @test nin(v3) == [nout(v1), nout(v2)] == [9, 11]

            @test_logs (:warn, r"Could not change nout of .* by -3! Relaxing constraints...") Δnout!(v3, -3)
            @test [nout(v1)] == nin(v2) == [9]
            @test nin(v3) == [nout(v1), nout(v2)] == [9, 9]
        end

    end

    @testset "Log size change" begin
        import Logging
        using Logging: Info
        using NaiveNASlib: NeuronIndices

        @testset "With strategy" begin
            av(in, outsize, name="av") = absorbvertex(SizeDummy(nout(in), outsize), in; traitdecoration=tf(name))

            v1 = inpt(3, "v1")
            v2 = av(v1, 10, "v2")
            v3 = av(v2, 4, "v3")

            Δninlog!(args...) = Δsize!(logafterΔsize(;level=Info, base=ΔNin(args...)))
            Δnoutlog!(args...) = Δsize!(logafterΔsize(;level=Info, base=ΔNout(args...)))

            Δnoutloginds!(args...) = Δsize!(NeuronIndices(), logafterΔsize(;level=Info, base=ΔNout(args...)))

            @test @test_logs (:info, "Change nin of v3 by 3") (:info, "Change nout of v2 by 3")  match_mode=:any Δninlog!(v3 => 3)

            @test @test_logs (:info, "Change nout of v2 by -3") (:info, "Change nin of v3 by -3") Δnoutlog!(v2 => -3)

            v4 = av(v1, nout(v3), "v4")
            v5 = "v5" >> v3 + v4

            @test @test_logs (:info, "Change nin of v5 by 30 and 30") (:info, "Change nout of v5 by 30") (:info, "Change nout of v3 by 30") (:info, "Change nout of v4 by 30") Δnoutlog!(v5 => 30)

            @test @test_logs (:info, "Change nin of v5 by [1,…, 34, -1×10] and [1,…, 34, -1×10]") (:info, "Change nout of v5 by [1,…, 34, -1×10]") (:info, "Change nout of v3 by [1,…, 34, -1×10]") (:info, "Change nout of v4 by [1,…, 34, -1×10]") Δnoutloginds!(v5 => 10)
        end

        @testset "With trait" begin
            traitfun(name) = t -> AfterΔSizeTrait(logafterΔsize(;level=Logging.Info), NamedTrait(name, t))
            avlogged(in, size ;name = "av") = absorbvertex(SizeDummy(nout(in), size), in; traitdecoration=traitfun(name))

            v1 = inpt(3, "v1")
            v2 = avlogged(v1, 10, name="v2")
            v3 = avlogged(v2, 4, name="v3")

            @test @test_logs (:info, "Change nin of v3 by 3") (:info, "Change nout of v2 by 3")  match_mode=:any Δnin!(v3, 3)

            @test @test_logs (:info, "Change nout of v2 by -3") (:info, "Change nin of v3 by -3") Δnout!(v2, -3)

            v4 = avlogged(v1, nout(v3), name="v4")
            v5 = traitconf(traitfun("v5")) >> v3 + v4

            @test @test_logs (:info, "Change nin of v5 by 30 and 30") (:info, "Change nout of v5 by 30") (:info, "Change nout of v3 by 30") (:info, "Change nout of v4 by 30") Δnout!(v5, 30)

            @test @test_logs (:info, "Change nin of v5 by [1,…, 34, -1×10] and [1,…, 34, -1×10]") (:info, "Change nout of v5 by [1,…, 34, -1×10]") (:info, "Change nout of v3 by [1,…, 34, -1×10]") (:info, "Change nout of v4 by [1,…, 34, -1×10]") Δsize!(NeuronIndices(), ΔNoutExact(v5, 10))
        end
    end

    @testset "WithKwargs" begin
        mutable struct WantKwargsSize
            insize::Int
            outsize::Int
        end
        NaiveNASlib.nout(w::WantKwargsSize) = w.outsize
        NaiveNASlib.nin(w::WantKwargsSize) = [w.insize]
        NaiveNASlib.Δsizetype(::WantKwargsSize) = NaiveNASlib.ScalarSize()

        seenkwargs = Dict()
        function NaiveNASlib.Δsize!(w::WantKwargsSize, ins::AbstractVector{<:Integer}, outs::Integer; kwargs...)
            if ins[] !== missing
                w.insize = ins[]
            end
            w.outsize =outs
            merge!(seenkwargs, kwargs)
        end

        iv = inpt(3)
        v1 = absorbvertex(WantKwargsSize(3, 5), iv; traitdecoration=tf("v1"))

        @test Δsize!(WithKwargs(ΔNout(v1 => 2); a=3, b=4)) == true

        @test nout(v1) == 7
        @test seenkwargs == Dict(:a => 3, :b => 4)

        @test Δsize!(validateafterΔsize(base=WithKwargs(ΔNout(v1 =>-2); c = 5))) == true
        @test nout(v1) == 5
        @test seenkwargs == Dict(:a => 3, :b => 4, :c=>5)
    end
end


text:
const data = parse.(Int, split(readline(), ","))

# Part One

minimum(UnitRange(extrema(data)...)) do i
    sum(abs.(data .- i))
end |> println

# Part Two

minimum(UnitRange(extrema(data)...)) do i
    sum(data) do d
        abs(d - i) * (abs(d - i) + 1) ÷ 2
    end
end |> println

ext: jl

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
import Term.style: apply_style
import Term: Panel

pts = """Lorem[red] ipsum dolor s[/red]it amet, consectetur adipiscing elit,
ed do eiusmod tempor [bold blue]incididu[underline]nt ut labore et dolore magna aliqua. Ut enim ad minim
 veniam, quis nos[/underline]trud exercitation ullamco laboris nisi ut aliquip ex
 ea commodo consequat. Duis aute [/bold blue][red on_black]irure dolor in reprehenderit
 in voluptate velit esse cillum dolore eu fugiat nulla
 pariatur. Excepteur sint occaecat[/red on_black] cupidatat non proident,
 sunt in [green]culpa qui officia[/green] deserunt mollit anim
 id est laborum."""

print(apply_style(pts))
@time apply_style(pts);

p2 = apply_style(pts)
@time print(Panel(pts; style="red"))


text:
using Optim

# Proxy function for optimization
function optFunc(x::Vector{CONST_TYPE}, dataset::Dataset{T}, baseline::T,
                 tree::Node, options::Options)::T where {T<:Real}
    setConstants(tree, x)
    return scoreFunc(dataset, baseline, tree, options)
end

# Use Nelder-Mead to optimize the constants in an equation
function optimizeConstants(dataset::Dataset{T},
                           baseline::T, member::PopMember,
                           options::Options)::PopMember where {T<:Real}

    nconst = countConstants(member.tree)
    if nconst == 0
        return member
    end
    x0 = getConstants(member.tree)
    f(x::Vector{CONST_TYPE})::T = optFunc(x, dataset, baseline, member.tree, options)
    if nconst == 1
        algorithm = Newton
    else
        algorithm = NelderMead
    end

    try
        result = optimize(f, x0, algorithm(), Optim.Options(iterations=100))
        # Try other initial conditions:
        for i=1:options.nrestarts
            new_start = x0 .* (convert(CONST_TYPE, 1) .+ convert(CONST_TYPE, 1//2)*randn(CONST_TYPE, size(x0)[1]))
            tmpresult = optimize(f, new_start, algorithm(), Optim.Options(iterations=100))

            if tmpresult.minimum < result.minimum
                result = tmpresult
            end
        end

        if Optim.converged(result)
            setConstants(member.tree, result.minimizer)
            member.score = convert(T, result.minimum)
            member.birth = getTime()
        else
            setConstants(member.tree, x0)
        end
    catch error
        # Fine if optimization encountered domain error, just return x0
        if isa(error, AssertionError)
            setConstants(member.tree, x0)
        else
            throw(error)
        end
    end
    return member
end


text:
import LinearAlgebra: ldiv!, rdiv!

ldiv!(α::Number, B::BlasArray) = lmul!(1 / α, B)
rdiv!(B::BlasArray, α::Number) = rmul!(B, 1 /α )

ldiv!(U::UniformScaling, B::BlasArray) = lmul!(1 / U.λ, B)
rdiv!(B::BlasArray, U::UniformScaling) = rmul!(B, 1 / U.λ)

# unwrap to-be-updated triangular numerators by sending ldiv! to !rdiv and vice-versa
# for (LoTr,UpTr) in ( (LowerTriangular, UpperTriangular),
#                      (UpperTriangular, LowerTriangular),
#                      (UnitLowerTriangular, UnitUpperTriangular),
#                      (UnitUpperTriangular, UnitLowerTriangular)
#                    )
#     for (tra,Tra) in ( (transpose,Transpose), (adjoint,Adjoint) )
#         @eval begin
#             ldiv!(L::$LoTr{T}, R::$Tra{T,<:$UpTr{T}}) where T =
#                 $tra(rdiv!($tra(R),$tra(L)))
#             #
#             rdiv!(L::$Tra{T,<:$LoTr{T}}, R::$UpTr{T}) where T =
#                 $tra(ldiv!($tra(R),$tra(L)))
#             #
#         end
#     end
# end


text:
export
    deltaT_seconds,
    dt_to_ut,
    tdt_to_bdt

# deltat_table is a list of tuples (jd, seconds), giving deltaT values
# for the beginnings of years in a historical range. [Meeus-1998: table 10.A]
deltat_table = Array[
    [cal_to_jd(1620, 1, 1), 121.0],
    [cal_to_jd(1622, 1, 1), 112.0],
    [cal_to_jd(1624, 1, 1), 103.0],
    [cal_to_jd(1626, 1, 1), 95.0],
    [cal_to_jd(1628, 1, 1), 88.0],

    [cal_to_jd(1630, 1, 1), 82.0],
    [cal_to_jd(1632, 1, 1), 77.0],
    [cal_to_jd(1634, 1, 1), 72.0],
    [cal_to_jd(1636, 1, 1), 68.0],
    [cal_to_jd(1638, 1, 1), 63.0],

    [cal_to_jd(1640, 1, 1), 60.0],
    [cal_to_jd(1642, 1, 1), 56.0],
    [cal_to_jd(1644, 1, 1), 53.0],
    [cal_to_jd(1646, 1, 1), 51.0],
    [cal_to_jd(1648, 1, 1), 48.0],

    [cal_to_jd(1650, 1, 1), 46.0],
    [cal_to_jd(1652, 1, 1), 44.0],
    [cal_to_jd(1654, 1, 1), 42.0],
    [cal_to_jd(1656, 1, 1), 40.0],
    [cal_to_jd(1658, 1, 1), 38.0],

    [cal_to_jd(1660, 1, 1), 35.0],
    [cal_to_jd(1662, 1, 1), 33.0],
    [cal_to_jd(1664, 1, 1), 31.0],
    [cal_to_jd(1666, 1, 1), 29.0],
    [cal_to_jd(1668, 1, 1), 26.0],

    [cal_to_jd(1670, 1, 1), 24.0],
    [cal_to_jd(1672, 1, 1), 22.0],
    [cal_to_jd(1674, 1, 1), 20.0],
    [cal_to_jd(1676, 1, 1), 28.0],
    [cal_to_jd(1678, 1, 1), 16.0],

    [cal_to_jd(1680, 1, 1), 14.0],
    [cal_to_jd(1682, 1, 1), 12.0],
    [cal_to_jd(1684, 1, 1), 11.0],
    [cal_to_jd(1686, 1, 1), 10.0],
    [cal_to_jd(1688, 1, 1), 9.0],

    [cal_to_jd(1690, 1, 1), 8.0],
    [cal_to_jd(1692, 1, 1), 7.0],
    [cal_to_jd(1694, 1, 1), 7.0],
    [cal_to_jd(1696, 1, 1), 7.0],
    [cal_to_jd(1698, 1, 1), 7.0],

    [cal_to_jd(1700, 1, 1), 7.0],
    [cal_to_jd(1702, 1, 1), 7.0],
    [cal_to_jd(1704, 1, 1), 8.0],
    [cal_to_jd(1706, 1, 1), 8.0],
    [cal_to_jd(1708, 1, 1), 9.0],

    [cal_to_jd(1710, 1, 1), 9.0],
    [cal_to_jd(1712, 1, 1), 9.0],
    [cal_to_jd(1714, 1, 1), 9.0],
    [cal_to_jd(1716, 1, 1), 9.0],
    [cal_to_jd(1718, 1, 1), 10.0],

    [cal_to_jd(1720, 1, 1), 10.0],
    [cal_to_jd(1722, 1, 1), 10.0],
    [cal_to_jd(1724, 1, 1), 10.0],
    [cal_to_jd(1726, 1, 1), 10.0],
    [cal_to_jd(1728, 1, 1), 10.0],

    [cal_to_jd(1730, 1, 1), 10.0],
    [cal_to_jd(1732, 1, 1), 10.0],
    [cal_to_jd(1734, 1, 1), 11.0],
    [cal_to_jd(1736, 1, 1), 11.0],
    [cal_to_jd(1738, 1, 1), 11.0],

    [cal_to_jd(1740, 1, 1), 11.0],
    [cal_to_jd(1742, 1, 1), 11.0],
    [cal_to_jd(1744, 1, 1), 12.0],
    [cal_to_jd(1746, 1, 1), 12.0],
    [cal_to_jd(1748, 1, 1), 12.0],

    [cal_to_jd(1750, 1, 1), 12.0],
    [cal_to_jd(1752, 1, 1), 13.0],
    [cal_to_jd(1754, 1, 1), 13.0],
    [cal_to_jd(1756, 1, 1), 13.0],
    [cal_to_jd(1758, 1, 1), 14.0],

    [cal_to_jd(1760, 1, 1), 14.0],
    [cal_to_jd(1762, 1, 1), 14.0],
    [cal_to_jd(1764, 1, 1), 14.0],
    [cal_to_jd(1766, 1, 1), 15.0],
    [cal_to_jd(1768, 1, 1), 15.0],

    [cal_to_jd(1770, 1, 1), 15.0],
    [cal_to_jd(1772, 1, 1), 15.0],
    [cal_to_jd(1774, 1, 1), 15.0],
    [cal_to_jd(1776, 1, 1), 16.0],
    [cal_to_jd(1778, 1, 1), 16.0],

    [cal_to_jd(1780, 1, 1), 16.0],
    [cal_to_jd(1782, 1, 1), 16.0],
    [cal_to_jd(1784, 1, 1), 16.0],
    [cal_to_jd(1786, 1, 1), 16.0],
    [cal_to_jd(1788, 1, 1), 16.0],

    [cal_to_jd(1790, 1, 1), 16.0],
    [cal_to_jd(1792, 1, 1), 15.0],
    [cal_to_jd(1794, 1, 1), 15.0],
    [cal_to_jd(1796, 1, 1), 14.0],
    [cal_to_jd(1798, 1, 1), 13.0],

    [cal_to_jd(1800, 1, 1), 13.1],
    [cal_to_jd(1802, 1, 1), 12.5],
    [cal_to_jd(1804, 1, 1), 12.2],
    [cal_to_jd(1806, 1, 1), 12.0],
    [cal_to_jd(1808, 1, 1), 12.0],

    [cal_to_jd(1810, 1, 1), 12.0],
    [cal_to_jd(1812, 1, 1), 12.0],
    [cal_to_jd(1814, 1, 1), 12.0],
    [cal_to_jd(1816, 1, 1), 12.0],
    [cal_to_jd(1818, 1, 1), 11.9],

    [cal_to_jd(1820, 1, 1), 11.6],
    [cal_to_jd(1822, 1, 1), 11.0],
    [cal_to_jd(1824, 1, 1), 10.2],
    [cal_to_jd(1826, 1, 1), 9.2],
    [cal_to_jd(1828, 1, 1), 8.2],

    [cal_to_jd(1830, 1, 1), 7.1],
    [cal_to_jd(1832, 1, 1), 6.2],
    [cal_to_jd(1834, 1, 1), 5.6],
    [cal_to_jd(1836, 1, 1), 5.4],
    [cal_to_jd(1838, 1, 1), 5.3],

    [cal_to_jd(1840, 1, 1), 5.4],
    [cal_to_jd(1842, 1, 1), 5.6],
    [cal_to_jd(1844, 1, 1), 5.9],
    [cal_to_jd(1846, 1, 1), 6.2],
    [cal_to_jd(1848, 1, 1), 6.5],

    [cal_to_jd(1850, 1, 1), 6.8],
    [cal_to_jd(1852, 1, 1), 7.1],
    [cal_to_jd(1854, 1, 1), 7.3],
    [cal_to_jd(1856, 1, 1), 7.5],
    [cal_to_jd(1858, 1, 1), 7.6],

    [cal_to_jd(1860, 1, 1), 7.7],
    [cal_to_jd(1862, 1, 1), 7.3],
    [cal_to_jd(1864, 1, 1), 6.2],
    [cal_to_jd(1866, 1, 1), 5.2],
    [cal_to_jd(1868, 1, 1), 2.7],

    [cal_to_jd(1870, 1, 1), 1.4],
    [cal_to_jd(1872, 1, 1), -1.2],
    [cal_to_jd(1874, 1, 1), -2.8],
    [cal_to_jd(1876, 1, 1), -3.8],
    [cal_to_jd(1878, 1, 1), -4.8],

    [cal_to_jd(1880, 1, 1), -5.5],
    [cal_to_jd(1882, 1, 1), -5.3],
    [cal_to_jd(1884, 1, 1), -5.6],
    [cal_to_jd(1886, 1, 1), -5.7],
    [cal_to_jd(1888, 1, 1), -5.9],

    [cal_to_jd(1890, 1, 1), -6.0],
    [cal_to_jd(1892, 1, 1), -6.3],
    [cal_to_jd(1894, 1, 1), -6.5],
    [cal_to_jd(1896, 1, 1), -6.2],
    [cal_to_jd(1898, 1, 1), -4.7],

    [cal_to_jd(1900, 1, 1), -2.8],
    [cal_to_jd(1902, 1, 1), -0.1],
    [cal_to_jd(1904, 1, 1), 2.6],
    [cal_to_jd(1906, 1, 1), 5.3],
    [cal_to_jd(1908, 1, 1), 7.7],

    [cal_to_jd(1910, 1, 1), 10.4],
    [cal_to_jd(1912, 1, 1), 13.3],
    [cal_to_jd(1914, 1, 1), 16.0],
    [cal_to_jd(1916, 1, 1), 18.2],
    [cal_to_jd(1918, 1, 1), 20.2],

    [cal_to_jd(1920, 1, 1), 21.1],
    [cal_to_jd(1922, 1, 1), 22.4],
    [cal_to_jd(1924, 1, 1), 23.5],
    [cal_to_jd(1926, 1, 1), 23.8],
    [cal_to_jd(1928, 1, 1), 24.3],

    [cal_to_jd(1930, 1, 1), 24.0],
    [cal_to_jd(1932, 1, 1), 23.9],
    [cal_to_jd(1934, 1, 1), 23.9],
    [cal_to_jd(1936, 1, 1), 23.7],
    [cal_to_jd(1938, 1, 1), 24.0],

    [cal_to_jd(1940, 1, 1), 24.3],
    [cal_to_jd(1942, 1, 1), 25.3],
    [cal_to_jd(1944, 1, 1), 26.2],
    [cal_to_jd(1946, 1, 1), 27.3],
    [cal_to_jd(1948, 1, 1), 28.2],

    [cal_to_jd(1950, 1, 1), 29.1],
    [cal_to_jd(1952, 1, 1), 30.0],
    [cal_to_jd(1954, 1, 1), 30.7],
    [cal_to_jd(1956, 1, 1), 31.4],
    [cal_to_jd(1958, 1, 1), 32.2],

    [cal_to_jd(1960, 1, 1), 33.1],
    [cal_to_jd(1962, 1, 1), 34.0],
    [cal_to_jd(1964, 1, 1), 35.0],
    [cal_to_jd(1966, 1, 1), 36.5],
    [cal_to_jd(1968, 1, 1), 38.3],

    [cal_to_jd(1970, 1, 1), 40.2],
    [cal_to_jd(1972, 1, 1), 42.2],
    [cal_to_jd(1974, 1, 1), 44.5],
    [cal_to_jd(1976, 1, 1), 46.5],
    [cal_to_jd(1978, 1, 1), 48.5],

    [cal_to_jd(1980, 1, 1), 50.5],
    [cal_to_jd(1982, 1, 1), 52.2],
    [cal_to_jd(1984, 1, 1), 53.8],
    [cal_to_jd(1986, 1, 1), 54.9],
    [cal_to_jd(1988, 1, 1), 55.8],

    [cal_to_jd(1990, 1, 1), 56.9],
    [cal_to_jd(1992, 1, 1), 58.3],
    [cal_to_jd(1994, 1, 1), 60.0],
    [cal_to_jd(1996, 1, 1), 61.6],
    [cal_to_jd(1998, 1, 1), 63.0],

    #  the following are not from Meeus, but are taken from
    # http://hpiers.obspm.fr/eop-pc/earthor/ut1lod/ut1-tai.html
    [cal_to_jd(2000, 1, 1), 32.184 + 31.6445],
    [cal_to_jd(2002, 1, 1), 32.184 + 32.1158],
    [cal_to_jd(2004, 1, 1), 32.184 + 32.3896],
    [cal_to_jd(2006, 1, 1), 32.184 + 32.6612],
    [cal_to_jd(2008, 1, 1), 32.184 + 33.2733],
    [cal_to_jd(2010, 1, 1), 32.184 + 33.897],
    [cal_to_jd(2012, 1, 1), 32.184 + 34.511]
]

deltat_table_start = 1620
deltat_table_end = 2012

"""
    deltaT_seconds(jd)

Return deltaT as seconds of time.

Functions which calculate the deltaT correction to convert between
dynamical and universal time.

Reference: Jean Meeus, _Astronomical Algorithms_, second edition, 1998

For a historical range from 1620 to a recent year, we interpolate from a
table of observed values. Outside that range we use formulae.

Parameters:
    jd : Julian Day number
Returns:
    deltaT in seconds
"""
function deltaT_seconds(jd)
    yr, mo, day = jd_to_cal(jd)
    jd1, jd0, secs1, secs0 = (0,0,0,0)
    #
    # 1620 - 20xx - do linear interpolation
    #
    if yr == deltat_table_end
    	return deltat_table_end[deltat_table][1]
    end

    if deltat_table_start <= yr && yr < deltat_table_end

    # linear search to find entries that bracket our target date. We could
    # improve this with a binary search or some sort of index.
    # TODO deltaT_seconds() improve crappy code
    for (i, j) in enumerate(deltat_table[1:end-1])
      if jd >= deltat_table[i][1] && jd < deltat_table[i+1][1]
        jd1, secs1 = deltat_table[i+1][1:2]
        jd0, secs0 = deltat_table[i][1:2]
      end
    end
    #  simple linear interpolation between two values
    return ((jd - jd0) * (secs1 - secs0) / (jd1 - jd0)) + secs0
    end

    t = (yr - 2000) / 100.0

    #
    # before 948 [Meeus-1998: equation 10.1]
    #
    if yr < 948
      return polynomial([2177, 497, 44.1], t)
    end

    #
    # 948 - 1620 and after 2012 [Meeus-1998: equation 10.2)
    #
    result = polynomial([102, 102, 25.3], t)

    #
    # correction for 2012-2100 [Meeus-1998: pg 78]
    #
    if deltat_table_end <= yr && yr < 2100
      result = result + 0.37 * (yr - 2100)
    end
    return result
end

"""
    dt_to_ut(jd)

Convert Julian Day from dynamical to terrestrial universal time.

Parameters:
    jd : Julian Day number (dynamical time)
Returns:
    Julian Day number (universal time)
"""
function dt_to_ut(jd)
    global seconds_per_day
    return jd - deltaT_seconds(jd) / seconds_per_day
end

"""
    tdt_to_bdt(tt, jd)

Converts from terrestrial dynamic time to barycentric dynamic time

Parameters:
    tt : terrestrial dynamic time
    jd : Julian Day number (terrestrial dynamic time)

Returns:
    barycentric dynamic time
"""
function tdt_to_bdt(tt, jd)
    g = 357.33 + 0.9856003*(jd-2451545)
    tdb = tt + 0.001658*sin(g) + 0.000014*sin(2*g)
    return tdb
end


text:
"""
    _serialized_change_unwrap(fn::Symbol)

Creates a simple wrapper structure that calls a function transparently on the
internal precached model. Internal type is returned (because this would break
the consistency of serialization).
"""
macro _serialized_change_unwrap(fn::Symbol)
    docstring = """
        $fn(model::Serialized, ...)

    Calls [`$fn`](@ref) of the internal serialized model type.
    Returns the modified un-serialized model.
    """
    :(@doc $docstring $fn(model::Serialized, args...; kwargs...) =
        $fn(unwrap_serialized(model), args...; kwargs...))
end


text:
include("advection.jl")
function timegradient(  m:: BurgerModel, x :: BurgerData ,TR :: TruneNone)
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(x.n))

  return BurgerData(x.n,diffusion + advectionLL)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneAllS )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))

  return BurgerData(x.n,diffusion + advectionLL)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneSS )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionSS = advection(m, x, ResolveSS(TR.n_off))

  return BurgerData(x.n,diffusion + advectionLL + advectionSS)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneSV )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionLS = advection(m, x, ResolveLS(TR.n_off))
  advectionSL = advection(m, x, ResolveSL(TR.n_off))

  return BurgerData(x.n,diffusion + advectionLL + advectionSL + advectionLS)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneExternalBurgerSS )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionSS = advection(m,x, ExternalResolveBurgerSS(TR.n_off, TR.x))

  return BurgerData(x.n,diffusion + advectionLL + advectionSS)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneExternalBurgerSV )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionLS = advection(m,x, ExternalResolveBurgerLS(TR.n_off, TR.x))
  advectionSL = advection(m,x, ExternalResolveBurgerSL(TR.n_off, TR.x))

  return BurgerData(x.n,diffusion + advectionLL + advectionSL + advectionLS)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneExternalBurgerSVSS )
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionSS = advection(m,x, ExternalResolveBurgerSS(TR.n_off, TR.x))
  advectionLS = advection(m,x, ExternalResolveBurgerLS(TR.n_off, TR.x))
  advectionSL = advection(m,x, ExternalResolveBurgerSL(TR.n_off, TR.x))
  return BurgerData(x.n,diffusion + advectionLL + advectionSL + advectionLS + advectionSS)
end

function timegradient( m :: BurgerModel, x:: BurgerData, TR :: TruneExternalBurgerSSDynamicSV)
  xx = x.x
  diffusion =  - m.d * ([1:x.n;]).^2 .* xx
  advectionLL = advection(m, x, ResolveLL(TR.n_off))
  advectionSS = advection(m,x, ExternalResolveBurgerSS(TR.n_off, TR.x))
  advectionLS = advection(m, x, ResolveLS(TR.n_off))
  advectionSL = advection(m, x, ResolveSL(TR.n_off))
  return BurgerData(x.n,diffusion + advectionLL + advectionSL + advectionLS + advectionSS)
end

export timegradient


text:
import Base.show
using Random

mutable struct HelpModel <: CmdStanModels
    @shared_fields_stanmodels
    method::Help
end

function HelpModel(
  name::AbstractString,
  model::AbstractString,
  n_chains=[4];
  seed = RandomSeed(),
  init = Init(),
  output = Output(),
  tmpdir = mktempdir(),
  method = Help(:sample),
  kwargs...)

  !isdir(tmpdir) && mkdir(tmpdir)

  StanBase.update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))

  output_base = joinpath(tmpdir, name)
  exec_path = executable_path(output_base)
  cmdstan_home = get_cmdstan_home()

  error_output = IOBuffer()
  is_ok = cd(cmdstan_home) do
      success(pipeline(`make -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`;
                       stderr = error_output))
  end
  if !is_ok
      throw(StanModelError(model, String(take!(error_output))))
  end

  HelpModel(name, model, n_chains, seed, init, output,
    tmpdir, output_base, exec_path, String[], String[],
    Cmd[], String[], String[], String[], false, false,
    cmdstan_home, method)
end

function help_model_show(io::IO, m, compact::Bool)
  println(io, "  name =                    \"$(m.name)\"")
  println(io, "  method =                  $(m.method)")
  println(io, "  n_chains =                $(get_n_chains(m))")
  println(io, "  tmpdir =                  \"$(m.tmpdir)\"")
end

show(io::IO, m::HelpModel) = help_model_show(io, m, false)


text:
using
    Oceananigans,
    Oceananigans.Operators,
    Oceananigans.Grids,
    Oceananigans.Solvers,
    Oceananigans.Diagnostics,
    Oceananigans.OutputWriters,
    Oceananigans.TurbulenceClosures,
    Test,
    Random,
    JLD2,
    Printf,
    Statistics,
    OffsetArrays,
    FFTW

@hascuda using CuArrays

using Statistics: mean
using LinearAlgebra: norm
using GPUifyLoops: @launch, @loop
using NCDatasets: Dataset

using Oceananigans: architecture, device, launch_config, datatuples, with_tracers,
                    Face, Cell, interiorparent, location, TracerFields, fill_halo_regions!

import Oceananigans: interior, datatuple

using Oceananigans.Solvers: PoissonSolver, PPN, PNN, solve_poisson_3d!

using Oceananigans.Diagnostics: run_diagnostic, velocity_div!

using Oceananigans.TimeSteppers: _compute_w_from_continuity!

using Oceananigans.TurbulenceClosures

using Oceananigans.TurbulenceClosures: ∂x_caa, ∂x_faa, ∂x²_caa, ∂x²_faa,
                                       ∂y_aca, ∂y_afa, ∂y²_aca, ∂y²_afa,
                                       ∂z_aac, ∂z_aaf, ∂z²_aac, ∂z²_aaf,
                                       ▶x_caa, ▶x_faa, ▶y_aca, ▶y_afa,
                                       ▶z_aac, ▶z_aaf

using Oceananigans.AbstractOperations

using Oceananigans.AbstractOperations: Computation, compute!

const seed = 420  # Random seed to use for all pseudorandom number generators.

datatuple(A) = NamedTuple{propertynames(A)}(Array(data(a)) for a in A)

function get_output_tuple(output, iter, tuplename)
    file = jldopen(output.filepath, "r")
    output_tuple = file["timeseries/$tuplename/$iter"]
    close(file)
    return output_tuple
end

float_types = (Float32, Float64)

archs = (CPU(),)
@hascuda archs = (CPU(), GPU())

closures = (
            :ConstantIsotropicDiffusivity,
            :ConstantAnisotropicDiffusivity,
            :AnisotropicBiharmonicDiffusivity,
            :SmagorinskyLilly,
            :BlasiusSmagorinsky,
            :RozemaAnisotropicMinimumDissipation,
            :VerstappenAnisotropicMinimumDissipation
           )

@testset "Oceananigans" begin
    include("test_grids.jl")
    include("test_fields.jl")
    include("test_halo_regions.jl")
    include("test_operators.jl")
    include("test_poisson_solvers.jl")
    include("test_coriolis.jl")
    include("test_surface_waves.jl")
    include("test_buoyancy.jl")
    include("test_models.jl")
    include("test_time_stepping.jl")
    include("test_boundary_conditions.jl")
    include("test_forcings.jl")
    include("test_turbulence_closures.jl")
    include("test_dynamics.jl")
    include("test_diagnostics.jl")
    include("test_output_writers.jl")
    include("test_regression.jl")
    include("test_examples.jl")
    include("test_abstract_operations.jl")
    include("test_verification.jl")
end


text:
# code for Video 4.5:
### RDatasets:

Pkg.add("RDatasets")
using RDatasets
data = dataset("datasets", "iris")

 RDatasets.packages()
33x2 DataFrames.DataFrame
| Row | Package        |
|-----|----------------|
| 1   | "COUNT"        |
| 2   | "Ecdat"        |
| 3   | "HSAUR"        |
| 4   | "HistData"     |
| 5   | "ISLR"         |
| 6   | "KMsurv"       |
| 7   | "MASS"         |
| 8   | "SASmixed"     |
| 9   | "Zelig"        |
...

julia> RDatasets.datasets("COUNT")
15x5 DataFrames.DataFrame
| Row | Package | Dataset      | Title        | Rows  | Columns |
|-----|---------|--------------|--------------|-------|---------|
| 1   | "COUNT" | "affairs"    | "affairs"    | 601   | 18      |
| 2   | "COUNT" | "azdrg112"   | "azdrg112"   | 1798  | 4       |
| 3   | "COUNT" | "azpro"      | "azpro"      | 3589  | 6       |
| 4   | "COUNT" | "badhealth"  | "badhealth"  | 1127  | 3       |
| 5   | "COUNT" | "fasttrakg"  | "fasttrakg"  | 15    | 9       |
| 6   | "COUNT" | "lbw"        | "lbw"        | 189   | 10      |
| 7   | "COUNT" | "lbwgrp"     | "lbwgrp"     | 6     | 7       |
| 8   | "COUNT" | "loomis"     | "loomis"     | 410   | 11      |
| 9   | "COUNT" | "mdvis"      | "mdvis"      | 2227  | 13      |
| 10  | "COUNT" | "medpar"     | "medpar"     | 1495  | 10      |
| 11  | "COUNT" | "rwm"        | "rwm"        | 27326 | 4       |
| 12  | "COUNT" | "rwm5yr"     | "rwm5yr"     | 19609 | 17      |
| 13  | "COUNT" | "ships"      | "ships"      | 40    | 7       |
| 14  | "COUNT" | "titanic"    | "titanic"    | 1316  | 4       |
| 15  | "COUNT" | "titanicgrp" | "titanicgrp" | 12    | 5       |

### reading in .rda files:
# R code to make iris.rda:
# iris = read.table("iris.csv", header=TRUE, sep = ",")
# save(iris, file = "iris.rda")

using DataArrays, DataFrames
rdata = read_rda("iris.rda")
Dict{ASCIIString,Any} with 1 entry:
  "iris" => DataFrames.RVector{Any,0x13}(Any[DataFrames.RVector{Float64,0x0e}([5…

### RCall:
Pkg.add("RCall")
using RCall

# rprint:
rprint(:iris)
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
4            4.6         3.1          1.5         0.2     setosa
5            5.0         3.6          1.4         0.2     setosa
6            5.4         3.9          1.7         0.4     setosa
...
148          6.5         3.0          5.2         2.0  virginica
149          6.2         3.4          5.4         2.3  virginica
150          5.9         3.0          5.1         1.8  virginica

# global environment:
g = globalEnv
RCall.RObject{RCall.EnvSxp}
<environment: R_GlobalEnv>

julia> globalEnv[:z] = [collect(1:10)]
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> rprint(:z)
 [1]  1  2  3  4  5  6  7  8  9 10

# R types:
RCall.CharSxp   # character string
RCall.StrSxp    # character string vector
RCall.IntSxp    # integer vector
RCall.LglSxp    # logical vector
RCall.RealSxp   # real vector
RCall.VecSxp    # list
RCall.CplxSxp   # complex vector
RCall.ClosSxp   # function closure
RCall.Sxp       # symbolic expression
RCall.NilSxp    # R's null value
RCall.RObject   # Julia wrapper for an R object

julia> RObject(10:20)
RCall.RObject{RCall.IntSxp}
 [1] 10 11 12 13 14 15 16 17 18 19 20

# reval:
julia> reval("0.1:1.0")
RCall.RObject{RCall.RealSxp}
[1] 0.1

julia> reval("'Julia is nice!'")
RCall.RObject{RCall.StrSxp}
[1] "Julia is nice!"

julia> reval("x <- list(val = 3.14)")
RCall.RObject{RCall.VecSxp}
$val
[1] 3.14

reval("plot(1:100)");
reval("dev.off()")

# rcall:
x = [0:0.01:pi]
y = cos(x)
rprint(rcall(:plot, x, y, xlab = "x", ylab = "cos(x)"))
reval("dev.off()")

# rcopy:
julia>rcopy("V <- runif(100, 1, 5)")
100-element Array{Float64,1}:
 3.06808
 4.6234
 4.5018
 ...
 julia>rcopy("hist(V)")
 julia>reval("dev.off()")
 RCall.RObject{RCall.IntSxp}
null device
          1

julia> rcopy("ls()")
4-element Array{ASCIIString,1}:
 "dfls"
 "V"
 "x"
 "z"

 julia> rprint(:x)
$val
[1] 3.14

julia> rprint(:z)
 [1]  1  2  3  4  5  6  7  8  9 10

julia> rprint(:V)
  [1] 3.068083 4.623404 4.501800 2.437669 4.128898 4.861698 3.012756 4.891440
  [9] 3.350561 4.791128 2.049340 1.947154 3.770557 1.221948 1.007942 3.359711
 [17] 4.338644 2.712553 4.263887 2.711554 2.945623 1.717035 2.083531 2.262899
 [25] 1.186694 4.646604 2.155006 2.757616 1.425041 1.394638 2.476713 3.235536
 [33] 4.356797 2.773285 2.885660 2.617869 4.587544 2.731121 4.124488 4.072608
 [41] 2.442061 2.058243 1.829554 3.372161 4.750352 2.410379 4.970044 3.017737
 [49] 2.910501 2.998271 1.672793 2.152949 2.292460 1.880926 4.081929 1.008254
 [57] 3.206837 4.826909 3.350008 1.040392 2.601423 1.249924 3.210415 2.971683
 [65] 2.951253 3.452599 3.583535 1.142748 3.778362 1.699750 3.909784 4.060263
 [73] 3.054700 2.055832 3.026260 3.359913 2.419350 4.900414 2.526388 1.627130
 [81] 4.716158 2.262563 3.905570 4.149724 3.466545 1.327037 3.255447 4.308987
 [89] 3.593277 1.124265 4.675755 1.463580 1.191946 2.303455 4.885928 3.767980
 [97] 2.055996 4.927590 4.352627 4.888622

julia> rprint(:dfls)
NULL

max_stars_repo_path: Section 4/4_5.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 4/4_5.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 4/4_5.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
# Functions for calculating zero variance (ZV) MCMC estimators, see
# Mira A, Solgi R, Imparato D. Zero Variance Markov Chain Monte Carlo for Bayesian Estimators. Statistics and
# Computing, 2013, 23 (5), pp 653-662

export linearZv, quadraticZv

# Functions for calculating ZV-MCMC estimators using linear polynomial
function linearZv(chain::Matrix{Float64}, grad::Matrix{Float64})
  npars = size(chain, 2)

  covAll = Array(Float64, npars+1, npars+1, npars)
  precision = Array(Float64, npars, npars, npars)
  sigma = Array(Float64, npars, npars)
  a = Array(Float64, npars, npars)

  z = -grad/2

  for i = 1:npars
     covAll[:, :, i] = cov([z chain[:, i]])
     precision[:, :, i] = inv(covAll[1:npars, 1:npars, i])
     sigma[:, i] = covAll[1:npars, npars+1, i]
     a[:, i] = -precision[:, :, i]*sigma[:, i]
  end

  chain = chain+z*a

  return chain, a
end

function linearZv(c::MCMCChain)
  if isa(c.task.sampler, ARS)
    indx = find(c.diagnostics["accept"])
  else
    indx = 1:size(c.samples, 1)
  end
  linearZv(c.samples[indx, :], c.gradients[indx, :])
end

# Functions for calculating ZV-MCMC estimators using quadratic polynomial
function quadraticZv(chain::Matrix{Float64}, grad::Matrix{Float64})
  nsamples, npars = size(chain)
  k = convert(Int, npars*(npars+3)/2)
  l = 2*npars+1

  zQuadratic = Array(Float64, nsamples, k)
  covAll = Array(Float64, k+1, k+1, npars)
  precision = Array(Float64, k, k, npars)
  sigma = Array(Float64, k, npars)
  a = Array(Float64, k, npars)

  z = -grad/2

  zQuadratic[:, 1:npars] = z
  zQuadratic[:, (npars+1):(2*npars)] = 2*z.*chain.-1
  for i = 1:(npars-1)
    for j = (i+1):npars
      zQuadratic[:, l] = chain[:, i].*z[:, j]+chain[:, j].*z[:, i]
      l += 1
    end
  end

  for i = 1:npars
    covAll[:, :, i] = cov([zQuadratic chain[:, i]]);
    precision[:, :, i] = inv(covAll[1:k, 1:k, i])
    sigma[:, i] = covAll[1:k, k+1, i]
    a[:, i] = -precision[:, :, i]*sigma[:, i]
  end

  zvChain = chain+zQuadratic*a;

  return zvChain, a
end

function quadraticZv(c::MCMCChain)
  if isa(c.task.sampler, ARS)
    indx = find(c.diagnostics["accept"])
  else
    indx = 1:size(c.samples, 1)
  end
  quadraticZv(c.samples[indx, :], c.gradients[indx, :])
end


text:
using ApproximateGPs
using Test

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


text:
using StatsBase
include("VRP_Heuristic.jl")
include("PDI_heuristic_resolution.jl")
include("BranchAndCut.jl")

function getAllInstancesPath(directory = "./PRP_instances")
    return readdir(directory, join = true)
end

function getAllInstancesPathFromType(type = "A")
    if type != "A" && type != "B"
        println("Choisir un type entre A et B")
        return []
    else
        paths = getAllInstancesPath()
        return filter(e -> first(split(last(split(e, "/")), "_")) == type, paths) #instances_paths
    end
end

function getAllInstancesPathFromTypeAndSize(type = "A", size = 14)
    allSizes = [14, 50, 100, 200]
    if !in(size, allSizes)
        print("Pas d'instance de taille $size")
        return []
    else
        instances_path_rightType = getAllInstancesPathFromType(type)

        return filter(e -> parse(Int, split(last(split(e, "/")), "_")[2]) == size, instances_path_rightType) #instances_path_rightType_rightSize
    end
end
function _try_get(f, model)
    try
        return f(model)
    catch
        return missing
    end
end
function _get_solution_dict(model)
    dict = Dict{String,Float64}()
    if has_values(model)
        for x in all_variables(model)
            variable_name = name(x)
            if !isempty(variable_name)
                dict[variable_name] = value(x)
            end
        end
    end
    return dict
end
function _get_constraint_dict(model)
    dict = Dict{String,Float64}()
    if has_duals(model)
        for (F, S) in list_of_constraint_types(model)
            for constraint in all_constraints(model, F, S)
                constraint_name = name(constraint)
                if !isempty(constraint_name)
                    dict[constraint_name] = dual(constraint)
                end
            end
        end
    end
    return dict
end
function getSelectedInstances(nbA14 = 20, nbA50 = 3, nbA100 = 0, nbB50 = 3, nbB100 = 0, nbB200 = 0)
    A14 = getAllInstancesPathFromTypeAndSize("A", 14)
    A50 = getAllInstancesPathFromTypeAndSize("A", 50)
    A100 = getAllInstancesPathFromTypeAndSize("A", 100)
    B50 = getAllInstancesPathFromTypeAndSize("B", 50)
    B100 = getAllInstancesPathFromTypeAndSize("B", 100)
    B200 = getAllInstancesPathFromTypeAndSize("B", 200)
    return [sample(A14, nbA14), sample(A50, nbA50), sample(A100, nbA100), sample(B50, nbB50), sample(B100, nbB100), sample(B200, nbB200)]
end

function callBranchAndCut(instances)
    for instance in instances
        #get the instance
        params, nodes, demands, costs = readPRP(instance)

        #Optimize
        model = BranchAndCutPDI(params, nodes, demands, costs)

        #create the directory
        timestamp = Dates.format(now(), "YYYYmmdd-HHMMSS")
        name = first(split(last(split(instance, "/")), "."))
        dir_name = joinpath(@__DIR__,"evaluations","Branch&Cut_$(name)_" * "$timestamp")
        @assert !ispath(dir_name) "Somebody else already created the directory"
        mkpath(dir_name)

        #save the circuits
        allCircuits = PDItoCircuits(model, params, nodes, demands, costs)
        saveMultiCircuits(params, allCircuits, dir_name * "/all_circuits.png")

        #save the model
        #write_to_file(model, dir_name * "/model_$(name).mps")
        towrite=["solver_name",solver_name(model),"\ntermination_status :",
        termination_status(model),"\nprimal_status :",
        primal_status(model),"\ndual_status :",
        dual_status(model),"\nraw_status :",
        raw_status(model),"\nresult_count :",
        result_count(model),"\nhas_values :",
        has_values(model),"\nhas_duals :",
        has_duals(model),"\nobjective_value :",
        _try_get(objective_value, model),"\nobjective_bound :",
        _try_get(objective_bound, model),"\ndual_objective_value :",
        _try_get(dual_objective_value, model),"\nSolutions :",
        _get_solution_dict(model),"\nConstraints :",
        _get_constraint_dict(model),"\nsolve_time :",
        _try_get(solve_time, model),"\nsimplex_iterations :",
        _try_get(simplex_iterations, model),"\nbarrier_iterations :",
        _try_get(barrier_iterations, model),"\nnode_count :",
        _try_get(node_count, model)]
        s=""
        for w in towrite
            s=s*string(w)*"\n"
        end
        file = open(dir_name *"/resume_exec_$(name).txt", "w")
	    write(file, s)
	    close(file)
    end
end

function callVRPHeuristic_and_Exact(instances)
    for instance in instances
        #get the instance
        params, nodes, demands, costs = readPRP(instance)

        #create the directory
        timestamp = Dates.format(now(), "YYYYmmdd-HHMMSS")
        name = first(split(last(split(instance, "/")), "."))
        dir_name = joinpath(@__DIR__,"evaluations","VRPHeuristic_and_Exact$(name)_" * "$timestamp")
        @assert !ispath(dir_name) "Somebody else already created the directory"
        mkpath(dir_name)

        #open the files
        file_heuristique = open(dir_name *"/circuits_de_clarkWright$(name).txt", "w")
        file_VRP = open(dir_name *"/circuits_de_VRP$(name).txt", "w")
        file_model=open(dir_name *"/all_models_de_VRP$(name).txt", "w")
	    write(file_heuristique, "t\trealisable\ttime\tcircuit\n")
        write(file_VRP, "t\trealisable\ttime\tcircuit\n")

        #get the circuits from metaheuristic and VRP for all time steps and save the model for all time steps
        for t in 1:params["l"]
            #save the heuristic's circuit
            start=now()
            circuits_temps_t_heuristique = mixMetaheuristic(clark_wright(params, nodes, demands, costs,t),params,costs)
            elapsed=now()-start
            isRealisable=true
            for circuit in circuits_temps_t_heuristique
                isRealisable=isRealisable&&length(circuit)<=params["k"]
                if !isRealisable
                    break
                end
            end
            write(file_heuristique, "$t\t$isRealisable\t$elapsed\t$(circuits_temps_t_heuristique)\n")

            #save the VRP's circuits
            model = createVRP_MTZ(params, nodes, demands, costs, t)
            optimize!(model)
            circuits_temps_t_VRP =vrpToCircuit(model, params, false, t)
            isRealisabl_VRPe=length(circuits_temps_t_VRP)<=params["k"]
            write(file_VRP, "$t\t$isRealisabl_VRPe\t$(_try_get(solve_time, model))\t$(circuits_temps_t_VRP)\n")

            #save the model
            towrite=["solver_name",solver_name(model),"\ntermination_status :",
            termination_status(model),"\nprimal_status :",
            primal_status(model),"\ndual_status :",
            dual_status(model),"\nraw_status :",
            raw_status(model),"\nresult_count :",
            result_count(model),"\nhas_values :",
            has_values(model),"\nhas_duals :",
            has_duals(model),"\nobjective_value :",
            _try_get(objective_value, model),"\nobjective_bound :",
            _try_get(objective_bound, model),"\ndual_objective_value :",
            _try_get(dual_objective_value, model),"\nSolutions :",
            _get_solution_dict(model),"\nConstraints :",
            _get_constraint_dict(model),"\nsolve_time :",
            _try_get(solve_time, model),"\nsimplex_iterations :",
            _try_get(simplex_iterations, model),"\nbarrier_iterations :",
            _try_get(barrier_iterations, model),"\nnode_count :",
            _try_get(node_count, model)]
            s="\n#$t\n"
            for w in towrite
                s=s*string(w)*"\n"
            end
            write(file_model,s)
        end

	    close(file_heuristique)
        close(file_VRP)
        close(file_model)

    end
end

function callPDIHeuristic(instances,nbMaxIte, resoudreVRPwithHeuristic)
    for instance in instances
        for instance in instances
            #get the instance
            lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial=initialisation_PDI_heuristique(instance)

            #Optimize
            model,circuits=PDI_heuristique(lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial, nbMaxIte, resoudreVRPwithHeuristic)

            #create the directory
            timestamp = Dates.format(now(), "YYYYmmdd-HHMMSS")
            name = first(split(last(split(instance, "/")), "."))
            withVRP=resoudreVRPwithHeuristic ? "" : "withVRP"
            dir_name = joinpath(@__DIR__,"evaluations","Branch&Cut_Heuristic_$(withVRP)_$(name)_" * "$timestamp")
            @assert !ispath(dir_name) "Somebody else already created the directory"
            mkpath(dir_name)

            #save the circuits
            circuits_list=[[] for i in 1:length(circuits)]
            for (i,circuit) in circuits
                circuits_list[i]=circuit
            end
            saveMultiCircuits(params, circuits_list, dir_name * "/all_circuits.png")

            #save the model
            #write_to_file(model, dir_name * "/model_$(name).mps")
            towrite=["solver_name",solver_name(model),"\ntermination_status :",
            termination_status(model),"\nprimal_status :",
            primal_status(model),"\ndual_status :",
            dual_status(model),"\nraw_status :",
            raw_status(model),"\nresult_count :",
            result_count(model),"\nhas_values :",
            has_values(model),"\nhas_duals :",
            has_duals(model),"\nobjective_value :",
            _try_get(objective_value, model),"\nobjective_bound :",
            _try_get(objective_bound, model),"\ndual_objective_value :",
            _try_get(dual_objective_value, model),"\nSolutions :",
            _get_solution_dict(model),"\nConstraints :",
            _get_constraint_dict(model),"\nsolve_time :",
            _try_get(solve_time, model),"\nsimplex_iterations :",
            _try_get(simplex_iterations, model),"\nbarrier_iterations :",
            _try_get(barrier_iterations, model),"\nnode_count :",
            _try_get(node_count, model)]
            s=""
            for w in towrite
                s=s*string(w)*"\n"
            end
            file = open(dir_name *"/resume_exec_$(name).txt", "w")
            write(file, s)
            close(file)
        end
    end
end

function evaluatePDI_BranchAndCut(nbA14 = 30, nbA50 = 30, nbA100 = 10, nbB50 = 0, nbB100 = 0, nbB200 = 0) #time limit de 3h pour chaque instance
    allInstances = getSelectedInstances(nbA14, nbA50, nbA100, nbB50, nbB100, nbB200)
    for instances in allInstances
        if length(instances)!=0
            callBranchAndCut(instances)
        end
    end
end

function evaluateVRP_Heuristique_and_Exact(nbA14 = 0, nbA50 = 1, nbA100 = 0, nbB50 = 0, nbB100 = 0, nbB200 = 0)
    allInstances = getSelectedInstances(nbA14, nbA50, nbA100, nbB50, nbB100, nbB200)
    for instances in allInstances
        callVRPHeuristic_and_Exact(instances)
    end
end

function evaluatePDI_heuristique(nbA14 = 10, nbA50 = 5, nbA100 = 2, nbB50 = 0, nbB100 = 0, nbB200 = 0, nbMaxIte=100, resoudreVRPwithHeuristic=false)
    allInstances = getSelectedInstances(nbA14, nbA50, nbA100, nbB50, nbB100, nbB200)
    for instances in allInstances
        callPDIHeuristic(instances,nbMaxIte, resoudreVRPwithHeuristic)
    end
end

#evaluatePDI_BranchAndCut(20,0,0,0,0,0)
evaluateVRP_Heuristique_and_Exact()
#evaluatePDI_heuristique()

#= TODO :
1 - Comparer en terme de vitesse et en fitness sur des instances de tailles différentes:
    - la résolution de VRP avec
        - la formulation MTZ
        - les différentes métaheuristiques
2 - Analyser les performances du LSP en terme de vitesse et en fitness sur des instances de tailles différentes
3 - Comparer en terme de vitesse et en fitness sur des instances de tailles différentes:
    -  PDI avec la résolution
        - heuristique
        -  exacte
Donc à chaque fois :
- enregistrer le temps d'execution, la fitness, les solutions
- Répéter les executions sur toutes les instances de tailles de différentes (arreter si trop long)
- sauvegarder les résultats sur un .txt afin de pouvoir les plotter en python
=#


text:
using ProjB
using Test

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


text:
## This simple (and naive) implementation is useful to solve the closed-shell linear expansion hartree-fock equation

module QP232

  include("./types.jl")

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

mutable struct NSurf <: AbstractBasis{2}
    order_u :: Int
    order_v :: Int
    knots_u :: Vector{Float64}
    knots_v :: Vector{Float64}
    weights :: Matrix{Float64}
end

function NSurf()
    NSurf(1, 1,
        [-1.0, -1.0, 1.0, 1.0],
        [-1.0, -1.0, 1.0, 1.0],
        ones(2, 2))
end

function length(basis::NSurf)
    nu = length(basis.knots_u) - basis.order_u - 1
    nv = length(basis.knots_v) - basis.order_v - 1
    return nu*nv
end

function size(basis::NSurf)
    return (2, length(basis))
end

function eval_basis!(basis::NSurf, N::Vector, xi::Vec{2})
    pu = basis.order_u
    pv = basis.order_v
    tu = basis.knots_u
    tv = basis.knots_v
    w = basis.weights
    nu = length(tu)-pu-1
    nv = length(tv)-pv-1
    u, v = xi
    n = 1
    for i=1:nu
        for j=1:nv
            N[n] = w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)
            n += 1
        end
    end
    N ./= sum(N)
    return N
end


text:
@testset "Trie for exponents vectors" begin

    for _ in 1:10
        d = rand([i + 10 for i in 1:10])
        t = ExpVectTrie(d)
        push!(t, [0 for _ in 1:d])
        vects = [rand([i for i in 1:10], d) for _ in 1:20]
        for v in vects
            push!(t, v)
        end
        for v in vects
            diff, best = get_max_below(t, v)
            @test diff == 0
            @test best == v
        end

        svects = Set(vects)
        push!(svects, [0 for _ in 1:d])

        for _ in 1:20
            v = rand([i for i in 1:10], d)
            diff, best = get_max_below(t, v)
            if v in svects
                @test diff == 0
                @test v == best
            else
                @test (best in svects)
                @test diff == sum(v .- best)
            end
        end
    end

end


text:
"""
    struct FinishedChannel{T} <: AbstractChannel{T}

Used to store tasks which have been completed, but not yet removed from the
underlying `DAG`. Taking from an empty `FinishedChannel` will block.
"""
struct FinishedChannel{T} <: AbstractChannel{T}
    data::Vector{T}
    cond_take::Threads.Condition
    function FinishedChannel{T}() where {T}
        lock      = Threads.ReentrantLock()
        cond_take = Threads.Condition(lock)
        data      = Vector{T}()
        new(data,cond_take)
    end
end

Base.lock(c::FinishedChannel)   = lock(c.cond_take)
Base.unlock(c::FinishedChannel) = unlock(c.cond_take)

Base.length(c::FinishedChannel) = length(c.data)

function Base.take!(c::FinishedChannel)
    lock(c)
    try
        while isempty(c.data)
            wait(c.cond_take)
        end
        @trace "length_finished $(length(c.data)) $(time_ns())"
        v = popfirst!(c.data)
        @trace "length_finished $(length(c.data)) $(time_ns())"
        return v
    finally
        unlock(c)
    end
end

function Base.put!(c::FinishedChannel{T},t::T) where {T}
    lock(c)
    try
        @trace "length_finished $(length(c.data)) $(time_ns())"
        push!(c.data, t)
        @trace "length_finished $(length(c.data)) $(time_ns())"
        notify(c.cond_take)
    finally
        unlock(c)
    end
    return t
end

# https://discourse.julialang.org/t/how-to-kill-thread/34236/8
"""
    struct Stop

Singleton type used to safely interrupt a task reading from an `AbstractChannel.
"""
struct Stop end

const Stoppable{T} = Union{T,Stop}

"""
    abstract type TaskGraphScheduler

Structures implementing a strategy to evaluate a [`DAG`](@ref).

Concrete subtypes are expected to contain a `dag::DAG` field for storing the
task graph, and a `finished::AbstractChannel` field to keep track of completed
tasks. The interface requires the following methods:

-`spawn(t,sch)`
-`schedule(t,sch)`

## See also: [`JuliaScheduler`](@ref), [`PriorityScheduler`](@ref), [`StaticScheduler`](@ref)
"""
abstract type TaskGraphScheduler end

"""
    const SCHEDULER::Ref{TaskGraphScheduler}

The active scheduler being used.
"""
const SCHEDULER = Ref{TaskGraphScheduler}()

"""
    setscheduler!(r)

Set the active scheduler to `r`.
"""
setscheduler!(ex) = (SCHEDULER[] = ex)

"""
    getscheduler(sch)

Return the active scheduler.
"""
getscheduler() = (SCHEDULER[])

"""
    with_scheduler(f,sch)

Run `f`, but push `DataFlowTask`s to the scheduler `dag` in `sch` instead of the
default `dag`.
"""
function with_scheduler(f,sch)
    old = getscheduler()
    setscheduler!(sch)
    res = f()
    setscheduler!(old)
    return res
end

spawn(tj::DataFlowTask) = spawn(tj,getscheduler())
Base.schedule(tj::DataFlowTask) = schedule(tj,getscheduler())

addnode!(sch::TaskGraphScheduler,tj,check=true)     = addnode!(sch.dag,tj,check)
remove_node!(sch::TaskGraphScheduler,tj) = remove_node!(sch.dag,tj)

"""
    sync([sch::TaskGraphScheduler])

Wait for all nodes in `sch` to be finished before continuining. If called with
no arguments, use  the current scheduler.
"""
function sync(sch::TaskGraphScheduler=getscheduler())
    dag = sch.dag
    isempty(dag) || wait(dag.cond_empty)
    return sch
end

graphplot(sch::TaskGraphScheduler=getscheduler()) = graphplot(sch.dag)

"""
    struct JuliaScheduler{T} <: TaskGraphScheduler{T}

Implement a simple scheduling strategy which consists of delegating the
[`DataFlowTask`](@ref)s to the native Julia scheduler for execution immediately
after the data dependencies have been analyzed using its `dag::DAG`. This is
the default scheduler used by [`DataFlowTasks`](@ref).

The main advantage of this strategy is its simplicity and composability. The
main disadvantage is that there is little control over how the underlying `Task`s
are executed by the Julia scheduler (e.g., no priorities can be passed).

Calling `JuliaScheduler(sz)` creates a new scheduler with an empty `DAG` of
maximum capacity `sz`.
"""
mutable struct JuliaScheduler{T} <: TaskGraphScheduler
    dag::DAG{T}
    finished::FinishedChannel{Stoppable{T}}
    dag_worker::Task
    function JuliaScheduler{T}(sz = typemax(Int)) where {T}
        if sz <= 0
            throw(ArgumentError("Scheduler buffer size must be a positive integer"))
        end
        dag            = DAG{T}(sz)
        finished       = FinishedChannel{Stoppable{T}}()
        sch            = new(dag,finished)
        start_dag_worker(sch)
        return sch
    end
end
JuliaScheduler(args...) = JuliaScheduler{DataFlowTask}(args...)

"""
    start_dag_worker(sch)

Start a forever-running task associated with `sch` which takes nodes from
`finished` and removes them from the `dag`. The task blocks if `finished` is
empty.
"""
function start_dag_worker(sch::JuliaScheduler=getscheduler())
    task = @async while true
        @assert Threads.threadid() == 1 # sanity check?
        t = take!(sch.finished)
        t == Stop() && break
        # remove task `t` from the dag
        remove_node!(sch,t)
    end
    sch.dag_worker = errormonitor(task)
    return sch.dag_worker
end

function stop_dag_worker(sch::JuliaScheduler=getscheduler())
    t = sch.dag_worker
    @assert istaskstarted(t)
    if istaskdone(t)
        @warn "DAG worker already stopped"
    elseif istaskfailed(t)
        @warn "DAG worker failed"
        return sch.dag_worker
    else # expected result, task is running
        isempty(sch.dag) || @warn "Stopping DAG worker of a non-empty graph"
        put!(sch.finished, Stop())
    end
    return sch.dag_worker
end

# interface methods for JuliaScheduler
function spawn(tj,::JuliaScheduler)
    tj.task.sticky = false
    schedule(tj.task)
    return tj
end

function Base.schedule(tj,::JuliaScheduler)
    schedule(tj.task)
    return tj
end

function Base.show(io::IO, sch::JuliaScheduler)
    dag = sch.dag
    n = num_nodes(dag)
    e = num_edges(dag)
    f = length(sch.finished)
    s1 = n==1 ? "" : "s"
    s2 = f==1 ? "" : "s"
    s3 = e==1 ? "" : "s"
    print(io, typeof(sch)," with $n active node$s1, $f finished node$s2, and $e edge$s3 (capacity of $(dag.sz_max[]) nodes)")
end

"""
    macro dspawn expr data mode

Create a [`DataFlowTask`](@ref) and schedule it to run on any available thread.
The `data` and `mode` arguments are passed to the `DataFlowTask` constructor,
and can be used to indicate how the code in `expr` accesses `data`. These fields
are used to automatically infer task dependencies.

## Examples:

```jldoctest
using DataFlowTasks
using DataFlowTasks: R,W,RW

A = rand(5)

# create a task which writes to A
t1 = @dspawn begin
    sleep(1)
    fill!(A,0)
    println("finished writing")
end (A,) (W,)

# create a task which reads from A
t2 = @dspawn begin
    println("I automatically wait for `t1` to finish")
    sum(A)
end (A,) (R,)

fetch(t2) # 0

# output

finished writing
I automatically wait for `t1` to finish
0.0
```

Note that in the example above `t2` waited for `t1` because it read a data field
that `t1` accessed in a writtable manner.
"""
macro dspawn(expr,data,mode,p=0)
    quote
        t = @dtask $(esc(expr)) $(esc(data)) $(esc(mode)) $(esc(p))
        spawn(t)
    end
end

"""
    macro dasync(expr,data,mode)

Like [`@dspawn`](@ref), but schedules the task to run on the current thread.
"""
macro dasync(expr,data,mode,p=0)
    quote
        t = @dtask $(esc(expr)) $(esc(data)) $(esc(mode)) $(esc(p))
        schedule(t)
    end
end


text:
#####################
# Logical circuits that are structured,
# meaning that each conjunction is associated with a vtree node.
#####################

"Root of the structure logical circuit node hierarchy"
abstract type StructLogicalΔNode{V<:VtreeNode} <: LogicalΔNode end

"A structured logical leaf node"
abstract type StructLogicalLeafNode{V} <: StructLogicalΔNode{V} end

"A structured logical inner node"
abstract type StructLogicalInnerNode{V} <: StructLogicalΔNode{V} end

"A structured logical literal leaf node, representing the positive or negative literal of its variable"
mutable struct StructLiteralNode{V} <: StructLogicalLeafNode{V}
    literal::Lit
    vtree::V
    data
    bit::Bool
    StructLiteralNode{V}(l,v::V) where V = new{V}(l, v, nothing, false)
end

"""
A structured logical constant leaf node, representing true or false.
These are the only structured nodes that don't have an associated vtree node (cf. SDD file format)
"""
abstract type StructConstantNode{V} <: StructLogicalInnerNode{V} end

mutable struct StructTrueNode{V} <: StructConstantNode{V}
    data
    bit::Bool
    StructTrueNode{V}() where V = new{V}(nothing, false)
end

mutable struct StructFalseNode{V} <: StructConstantNode{V}
    data
    bit::Bool
    StructFalseNode{V}() where V = new{V}(nothing, false)
end

"A structured logical conjunction node"
mutable struct Struct⋀Node{V} <: StructLogicalInnerNode{V}
    children::Vector{<:StructLogicalΔNode{<:V}}
    vtree::V
    data
    bit::Bool
    Struct⋀Node{V}(c,v::V) where V = new{V}(convert(Vector{StructLogicalΔNode{V}}, c), v, nothing, false)
end

"A structured logical disjunction node"
mutable struct Struct⋁Node{V} <: StructLogicalInnerNode{V}
    children::Vector{<:StructLogicalΔNode{<:V}}
    vtree::V # could be leaf or inner
    data
    bit::Bool
    Struct⋁Node{V}(c,v::V) where V = new{V}(convert(Vector{StructLogicalΔNode{V}}, c), v, nothing, false)
end

HasVtree = Union{Struct⋁Node,Struct⋀Node,StructLiteralNode}

"A structured logical circuit represented as a bottom-up linear order of nodes"
const StructLogicalΔ{V} = AbstractVector{<:StructLogicalΔNode{V}}

#####################
# traits
#####################

@inline GateType(::Type{<:StructLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:StructConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Struct⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Struct⋁Node}) = ⋁Gate()

@inline node_type(::Type{<:StructLogicalΔNode}) = StructLogicalΔNode
#####################
# methods
#####################

@inline literal(n::StructLiteralNode)::Lit = n.literal
@inline constant(n::StructTrueNode)::Bool = true
@inline constant(n::StructFalseNode)::Bool = false
@inline children(n::StructLogicalInnerNode) = n.children

@inline vtree(n::HasVtree) = n.vtree

"Conjoin nodes in the same way as the example"
@inline function conjoin_like(example::StructLogicalΔNode, arguments::Vector)
    if isempty(arguments)
        StructTrueNode{PlainVtreeNode}()
    elseif example isa Struct⋀Node && children(example) == arguments
        example
    else
        Struct⋀Node{PlainVtreeNode}(arguments, vtree(example))
    end
end

"Disjoin nodes in the same way as the example"
@inline function disjoin_like(example::StructLogicalΔNode, arguments::Vector)
    if isempty(arguments)
        StructFalseNode{PlainVtreeNode}()
    elseif example isa Struct⋁Node && children(example) == arguments
        example
    else
        Struct⋁Node{PlainVtreeNode}(arguments, vtree(example))
    end
end

"Construct a new literal node like the given node's type"
literal_like(example::StructLogicalΔNode, lit::Lit) = StructLiteralNode{PlainVtreeNode}(lit, vtree(example))

@inline copy_node(n::Struct⋁Node, cns) = Struct⋁Node{PlainVtreeNode}(cns, vtree(n))
@inline copy_node(n::Struct⋀Node, cns) = Struct⋀Node{PlainVtreeNode}(cns, vtree(n))


text:
using JLD2
using Statistics

mutable struct TwoDGridWorldError <: LoggerKeyData
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int

    function TwoDGridWorldError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TwoDGridWorldSet.jld2") eval_set
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1

        new(zeros(4,num_logged_steps), eval_set, logger_init_info[LoggerInitKey.INTERVAL])
    end
end

function lg_start!(self::TwoDGridWorldError, env, agent)
    Q_est = hcat([get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)
    true_values = TwoDGridWorldError.get_true_values(env, self.eval_set["ests"])
    err = mean((Q_est - true_values) .^ 2, dims=2)
    self.error[:,1] = sqrt.(err)
end

function lg_step!(self::TwoDGridWorldError, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if (rem(cur_step_total, self.log_interval) == 0)
        ind = fld(cur_step_total, self.log_interval) + 1
        Q_est = hcat([get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)

        true_values = TwoDGridWorldError.get_true_values(env, self.eval_set["ests"])
        err = mean((Q_est - true_values) .^ 2, dims=2)

        self.error[:,ind] = sqrt.(err)
    end
end

function lg_episode_end!(self::TwoDGridWorldError, cur_step_in_episode, cur_step_total)
end

function save_log(self::TwoDGridWorldError, save_dict::Dict)
    save_dict[:twod_gridworld_error] = self.error
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters

using BenchmarkTools
using Profile
using Distributed
using CountPlus
using LinearAlgebra.BLAS: set_num_threads
using Base.Threads: nthreads

using ShaleDrillingLikelihood: test_parallel_simloglik!

num_i = 150
M = 100
nz = 51
nψ = 51

data_theta_sm, data_theta_lg = MakeTestData(; num_i=num_i, nz=nz, nψ=nψ)

nlogic_cores = length(Sys.cpu_info())
nphys_cores = Int(nlogic_cores/2)

rmprocs(filter(i -> i != 1, workers()))
nth = nthreads()
nproc = Int(nlogic_cores / nth)
pids = addprocs(nproc; enable_threaded_blas=true)

println("----------------------------------------------")
println("Using $nth threads and adding $nproc workers")
println("----------------------------------------------")

@everywhere begin
    using LinearAlgebra.BLAS: set_num_threads
    using Base.Threads: nthreads
    using ShaleDrillingLikelihood
    set_num_threads(nthreads())
end

data, theta = data_theta_sm

theta_peturb = 0.9 .* theta

leo = LocalEstObj(data,theta_peturb)
s = SimulationDraws(M, ShaleDrillingLikelihood.data(leo))
reo = RemoteEstObj(leo, M)
ew = EstimationWrapper(leo, reo)

stopcount!()
for dograd in false:true
    @eval @everywhere set_g_RemoteEstObj($reo)
    @btime test_parallel_simloglik!($ew, $theta_peturb, $dograd)
end

# end


text:
@testset "DataValueArray: Map" begin

using Test
using DataValues

# create m random arrays each with N dimensions and length dims[i] along
# dimension i for i=1:N
N = rand(2:5)
dims = Int[ rand(3:8) for i in 1:N ]
m = rand(4:6)
n = rand(1:3)
As = Array[ [rand(dims...) for i in 1:n] ; [rand(Int, dims...) for i in n+1:m] ]

Ms = [ rand(Bool, dims...) for i in 1:m ]
Xs = Array{DataValueArray, 1}()
for i in 1:m
    push!(Xs, DataValueArray(As[i]))
end
Ys = Array{DataValueArray, 1}()
for i in 1:m
    push!(Ys, DataValueArray(As[i], Ms[i]))
end

C = Array{Float64}(dims...)
Z = DataValueArray{Float64}(dims...)

R = map(|, Ms...)

f(x...) = sum(x)
dests = (C, Z)

# 1 arg
for (args, masks) in (
    ((As, Xs), fill((), m)), ((As, Ys), [ (Ms[i],) for i in 1:m ])
)
    for i in 1:m
        # map!
        map!(f, args[1][i], args[1][i]) # map!(f, As[i], As[i])
        map!(f, args[2][i], args[2][i]) # map!(f, Xs[i], Xs[i])
        @test isequal(args[2][i], DataValueArray(args[1][i], masks[i]...))
        # map
        A = map(f, args[1][i])
        X = map(f, args[2][i])
        @test isequal(X, DataValueArray(A, masks[i]...))
    end
end
# 2 arg
i, j = rand(1:m), rand(1:m)
S = map(|, Ms[i], Ms[j])
for (args, mask) in (
    ((As, Xs), ()), ((As, Ys), (S,))
)
    # map!
    map!(f, dests[1], args[1][i], args[1][j])
    map!(f, dests[2], args[2][i], args[2][j])
    @test isequal(dests[2], DataValueArray(dests[1], mask...))
    # map
    map(f, args[1][i], args[1][j])
    map(f, args[2][i], args[2][j])
    @test isequal(dests[2], DataValueArray(dests[1], mask...))
end
# n arg
for (args, mask) in (
    ((As, Xs), ()), ((As, Ys), (R,))
)
    # map!
    map!(f, dests[1], args[1]...)
    map!(f, dests[2], args[2]...)
    @test isequal(dests[2], DataValueArray(dests[1], mask...))
    # map
    map(f, args[1]...)
    map(f, args[2]...)
    @test isequal(dests[2], DataValueArray(dests[1], mask...))
end

# test map over empty DataValueArrays
X = DataValueArray(Int[])
h1(x) = 5.0*x
h2(x) = x
h2(x...) = prod(x)

Z1 = map(h1, X)
Z2 = map(h2, X)
@test isempty(Z1)
@test isa(Z1, DataValueArray{Float64})
@test isempty(Z2)
@test isa(Z2, DataValueArray{Int})

# if a function has no method for inner eltype of empty DataValueArray,
# result should be empty DataValueArray{Any}() for consistency with generic map()
h3(x::Float64...) = prod(x)
Z3 = map(h3, X)
@test isempty(Z3)
@test isa(Z3, DataValueArray{Any})
Z3 = map(h3, X, X)
@test isempty(Z3)
@test isa(Z3, DataValueVector{Any})
Z3 = map(h3, X, X, X)
@test isempty(Z3)
@test isa(Z3, DataValueArray{Any})

# test map over all null DataValueArray
n = rand(10:100)
Ys = [ DataValueArray(rand(Int, n), fill(true, n)) for i in 1:rand(3:5) ]

Z2 = map(h2, Ys[1])
@test isequal(Z2, DataValueArray{Int}(n))
@test isa(Z2, DataValueArray{Int})
Z2 = map(h2, Ys[1], Ys[2])
@test isequal(Z2, DataValueArray{Int}(n))
@test isa(Z2, DataValueArray{Int})
Z2 = map(h2, Ys...)
@test isequal(Z2, DataValueArray{Int}(n))
@test isa(Z2, DataValueArray{Int})

end