text:
@testset "Structure" begin
    import ONNXNaiveNASflux: OnnxNode, input, output, optype, array

    function serdeser(p::T, convfun = cfun(p)) where T
        iob = PipeBuffer();
        ONNX.writeproto(iob, p)
        return convfun(ONNX.readproto(iob, T()))
    end

    function serdeser(p::ONNX.ModelProto)
        iob = PipeBuffer();
        save(iob, p)
        return ONNXNaiveNASflux.extract(iob)
    end

    cfun(::ONNX.NodeProto) = np -> OnnxNode(np, ONNX.TensorProto[])
    cfun(::ONNX.TensorProto) = array
    cfun(::ONNX.GraphProto) = identity

    include("onnxruntime.jl")

    @testset "Nodes" begin
        using Statistics
        import ONNXNaiveNASflux: optype, actfuns, fluxlayers, invariantops
        using ONNXNaiveNASflux.NaiveNASflux
        import ONNXNaiveNASflux.NaiveNASflux: weights, bias
        import ONNXNaiveNASflux: AbstractProbe, nextname, newfrom, add!, genname, shape, nextshape
        struct NodeProbe{F, S} <: AbstractProbe
            name::String
            namefun::F
            shape::S
            protos::Vector{Any}
        end
        NodeProbe(name, namefun, shape=missing) = NodeProbe(name, namefun, shape, [])
        ONNXNaiveNASflux.add!(p::NodeProbe, n) = push!(p.protos, n)
        ONNXNaiveNASflux.nextname(p::NodeProbe) = p.namefun
        ONNXNaiveNASflux.newfrom(p::NodeProbe, pname, Δshape=identity) = NodeProbe(pname, p.namefun, nextshape(p, Δshape), p.protos)
        ONNXNaiveNASflux.newnamestrat(p::NodeProbe, f, pname = name(p)) = NodeProbe(pname, f, p.shape, p.protos)
        ONNXNaiveNASflux.name(p::NodeProbe) = p.name
        ONNXNaiveNASflux.shape(p::NodeProbe) = p.shape

        @testset "Paramfree op $(tc.op) attrs: $(pairs(tc.attr))" for tc in (
            (op=:Relu, attr = Dict(), fd=actfuns),
            (op=:Elu, attr = Dict(), fd=actfuns),
            (op=:Elu, attr = Dict(:alpha => 0.5f0), fd=actfuns),
            (op=:Selu, attr = Dict(), fd=actfuns),
            (op=:Selu, attr = Dict(:alpha => 1.5f0), fd=actfuns),
            (op=:Tanh, attr = Dict(), fd=actfuns),
            (op=:GlobalAveragePool, attr=Dict(), fd=invariantops),
            (op=:MaxPool, attr=Dict(:kernel_shape=>(1,2), :pads=>(2,1), :strides=>(2,2)), fd=fluxlayers),
            (op=:AveragePool, attr=Dict(:kernel_shape=>(3,2), :pads=>(1,0), :strides=>(2,2)), fd=fluxlayers),
            (op=:Dropout, attr=Dict(:ratio => 0.2f0), fd=fluxlayers),
            )

            inprobe = NodeProbe("input", f -> "output")

            outprobe = tc.fd[tc.op](tc.attr)(inprobe)

            @test length(outprobe.protos) == 1

            res = serdeser(outprobe.protos[])

            @test input(res) == [name(inprobe)]
            @test output(res) == [name(outprobe)]
            @test optype(res) == tc.op
            @test name(res) == name(outprobe)

            mexprev(v, x) = x
            mexprev(v, x::Tuple) = reverse(x)
            mexprev(::Val{:pads}, x::Tuple) = ONNXNaiveNASflux.padexpand(Val(length(x)), x)
            for (k,v) in tc.attr
                for (exp, act) in zip(mexprev(Val(k), v), res.attribute[k])
                    @test exp == act
                end
            end
        end

        @testset "Dims method $(tc.ot)" for tc in (
            (f=cat, dims=1, ndims=2, ot=:Concat, axname=:axis),
            (f=mean, dims=(2, 3), ndims=4, ot=:ReduceMean, axname=:axes),
            (f=dropdims, dims=(3,), ndims=3, ot=:Squeeze, axname=:axes)
            )
            inprobe = NodeProbe("input", f -> "output", Tuple(1:tc.ndims))

            outprobe = tc.f(inprobe, dims=tc.dims)

            @test length(outprobe.protos) == 1

            res = serdeser(outprobe.protos[])

            @test input(res) == [name(inprobe)]
            @test output(res) == [name(outprobe)]
            @test optype(res) == tc.ot
            @test name(res) == name(outprobe)
            expdims = tc.dims isa Tuple ? collect(tc.dims) : tc.dims
            @test ONNXNaiveNASflux.numpy2fluxdim.(res.attribute[tc.axname], tc.ndims) == expdims
        end

        @testset "Reshape" begin
            inprobe = NodeProbe("input", f -> "output", (:A, missing, 12))

            outprobe = reshape(inprobe, (0, 3, 2, Colon()))
            shapeout = shape(outprobe)

            @test length(outprobe.protos) == 2

            res = serdeser(outprobe.protos[1])
            newshape = serdeser(outprobe.protos[2])

            @test newshape == [-1, 2, 3, 0]

            @test input(res) == [name(inprobe), outprobe.protos[2].name]
            @test output(res) == [name(outprobe)]
            @test optype(res) == :Reshape
            @test name(res) == name(outprobe)
            @test ismissing(shapeout[end])
            @test collect(skipmissing(shapeout)) == [:A, 3, 2]
        end

        @testset "Flatten" begin
            import ONNXNaiveNASflux.NaiveNASflux.Flux: flatten
            inprobe = NodeProbe("input", f -> "output", (2, 3, 5))

            outprobe = Flux.flatten(inprobe)
            shapeout = shape(outprobe)

            @test length(outprobe.protos) == 1

            res = serdeser(outprobe.protos[1])

            @test input(res) == [name(inprobe)]
            @test output(res) == [name(outprobe)]
            @test optype(res) == :Flatten
            @test name(res) == name(outprobe)
            @test res.attribute[:axis] == -2

            op = ONNXNaiveNASflux.pseudotransparentops[optype(res)](res.attribute)
            indata = reshape(collect(1:2*3*5), 2,3,5)
            @test op(indata) == flatten(indata)
        end

        @testset "Pad expand" begin
            import ONNXNaiveNASflux: padexpand
            @test padexpand(Val(1), (1,)) == [1,1]
            @test padexpand(Val(2), (1,2)) == [2,1,2,1]
            @test padexpand(Val(3), (1,2,3)) == [3,2,1,3,2,1]

            @test padexpand(Val(1), (1,2)) == [1,2]
            @test padexpand(Val(2), (1,2,3,4)) == [3,1,4,2]
            @test padexpand(Val(3), (1,2,3,4,5,6)) == [5,3,1,6,4,2]
        end

        @testset "Layer with activation function $actfun" for actfun in (
            relu,
            elu,
            selu,
            tanh,
        )
            @testset "$(tc.layer) node" for tc in (
                (layer=Dense(3,4, actfun), indata=reshape(collect(Float32, 1:12), :, 4) .- 3),
                (layer=Conv((1,2), 3=>4, actfun; pad=(2,1), stride=(1,2), dilation=3), indata=reshape(collect(Float32, 1:2*3*9*9), 9,9,3,2) .- 5),
                (layer=Conv((2,3), 3=>4, actfun; pad=(1,2,3,4), stride=(1,2), dilation=3), indata=reshape(collect(Float32, 1:2*3*9*9), 9,9,3,2) .- 10),
                )

                inprobe = NodeProbe("input", genname, shape(layertype(tc.layer), nin(tc.layer)))

                outprobe = tc.layer(inprobe)

                @test length(outprobe.protos) == 4

                wp,bp,lp,ap = Tuple(outprobe.protos)

                ln = serdeser(lp)
                an = serdeser(ap)
                w = serdeser(wp)
                b = serdeser(bp)

                @test size(w) == size(weights(tc.layer))
                @test size(b) == size(bias(tc.layer))

                @test w ≈ ONNXNaiveNASflux.flipweights(layertype(tc.layer), weights(tc.layer))
                @test b ≈ bias(tc.layer)

                ln.attribute[:activation] = actfuns[Symbol(optype(an))](an.attribute)
                res = fluxlayers[optype(ln)](ln.attribute, w, b)

                resout = res(tc.indata)
                expout = tc.layer(tc.indata)

                @test size(resout) == size(expout)
                @test resout ≈ expout

                ortout, = onnxruntime_infer(tc.layer, tc.indata)
                @test size(ortout) == size(expout)
                @test ortout ≈ expout
            end
        end

        @testset "$(tc.layer) node no bias no act" for tc in (
            (layer=Dense(randn(Float32, 2,3), Flux.Zeros()), indata=reshape(collect(Float32, 1:12), :, 4) .- 3),
            (layer=Conv((1,1), 2=>3; bias=Flux.Zeros(3)), indata=reshape(collect(Float32, 1:2*3), 1,1,2,3) .- 3),
            )

            inprobe = NodeProbe("input", genname, shape(layertype(tc.layer), nin(tc.layer)))

            outprobe = tc.layer(inprobe)

            @test length(outprobe.protos) == 2

            wp,lp= Tuple(outprobe.protos)

            ln = serdeser(lp)
            w = serdeser(wp)

            @test size(w) == size(weights(tc.layer))
            @test w ≈ ONNXNaiveNASflux.flipweights(layertype(tc.layer), weights(tc.layer))

            res = fluxlayers[optype(ln)](ln.attribute, w)

            resout = res(tc.indata)
            expout = tc.layer(tc.indata)

            @test size(resout) == size(expout)
            @test resout ≈ expout

            ortout, = onnxruntime_infer(tc.layer, tc.indata)
            @test size(ortout) == size(expout)
            @test ortout ≈ expout
        end

        @testset "$(tc.layer) node" for tc in (
            (layer=RNN(3, 5, x -> Flux.elu(x, 0.1f0)), indata = reshape(collect(Float32, 1:12), :, 4) .- 3),
            (layer=LSTM(4, 3), indata = reshape(collect(Float32, 1:12), 4, :) .- 3),
            )
            import ONNXNaiveNASflux.NaiveNASflux: hiddenweights

            inprobe = NodeProbe("input", genname, shape(layertype(tc.layer), nin(tc.layer)))

            outprobe = tc.layer(inprobe)

            @test length(outprobe.protos) == 5

            wip,whp,bp,lp = Tuple(outprobe.protos)

            ln = serdeser(lp)
            wi = serdeser(wip)
            wh = serdeser(whp)
            b = serdeser(bp)

            res = fluxlayers[optype(ln)](ln.attribute, wi, wh, b)

            lt = layertype(tc.layer)
            @test size(weights(res)) == size(weights(tc.layer))
            @test size(hiddenweights(res)) == size(hiddenweights(tc.layer))
            @test size(bias(res)) == size(bias(res))

            @test weights(res) ≈ weights(tc.layer)
            @test hiddenweights(res) ≈ hiddenweights(tc.layer)
            @test bias(res) ≈ bias(res)

            resout = res(tc.indata)
            expout = tc.layer(tc.indata)

            @test size(resout) == size(expout)
            @test resout ≈ expout

            ortout, = onnxruntime_infer(tc.layer, reshape(tc.indata,size(tc.indata)...,1))
            ortout = dropdims(ortout; dims=3)

            @test size(ortout) == size(expout)
            @test ortout ≈ expout
        end

        @testset "$(tc.layer) node" for tc in (
            (layer=BatchNorm(3, relu; initβ = i -> collect(Float32, 1:i), initγ = i -> collect(Float32, i:-1:1), ϵ=1e-3, momentum = 0.78), indata=reshape(collect(Float32, 1:2*3*3), 2,3,3,1) .- 10),
            )

            inprobe = NodeProbe("input", genname)
            outprobe = tc.layer(inprobe)
            @test length(outprobe.protos) == 6

            ln, γ, β, μ, σ², an = Tuple(serdeser.(outprobe.protos))

            @test size(β) == size(tc.layer.β)
            @test size(γ) == size(tc.layer.γ)
            @test size(μ) == size(tc.layer.μ)
            @test size(σ²) == size(tc.layer.σ²)

            @test β ≈ tc.layer.β
            @test γ ≈ tc.layer.γ
            @test μ ≈ tc.layer.μ
            @test σ² ≈ tc.layer.σ²

            ln.attribute[:activation] = actfuns[Symbol(optype(an))](an.attribute)
            res = fluxlayers[optype(ln)](ln.attribute, γ, β, μ, σ²)

            @test string(res) == string(tc.layer)

            resout = res(tc.indata)
            expout = tc.layer(tc.indata)

            @test size(resout) == size(expout)
            @test resout ≈ expout

            ortout, = onnxruntime_infer(tc.layer, tc.indata)
            @test size(ortout) == size(expout)
            @test ortout ≈ expout
        end
    end

    @testset "Graphs" begin
        using ONNXNaiveNASflux.NaiveNASflux
        using ONNXNaiveNASflux: graphproto, modelproto, validate, AbstractVertex
        using NaiveNASflux: AbstractMutableComp

        dense(name, inpt::AbstractVertex, outsize, actfun=identity) = fluxvertex(name, Dense(nout(inpt), outsize, actfun), inpt)
        dense(inpt::AbstractVertex, outsize, actfun=identity) = fluxvertex(Dense(nout(inpt), outsize, actfun), inpt)

        convvertex(name, inpt::AbstractVertex, outsize, actfun=identity) = fluxvertex(name, Conv((1,1), nout(inpt) => outsize, actfun), inpt)

        bnvertex(name, inpt::AbstractVertex, actfun=identity) = fluxvertex(name, BatchNorm(nout(inpt), actfun), inpt)

        mpvertex(name, inpt::AbstractVertex) = fluxvertex(name, MaxPool((2,2); pad=(1,0), stride=(1,2)), inpt)

        fvertex(name, inpt::AbstractVertex, f) = invariantvertex(name, f, inpt)

        test_outputs(res::Tuple, exp, sizediff=0) = test_outputs(res[1], exp, sizediff)
        test_outputs(res::Tuple, exp::Tuple, sizediff=0) = test_outputs.(res, exp, sizediff)
        function test_outputs(res, exp, sizediff=0)
            # Fix for differences in recurrent shapes. Flux ises 2D and onnx uses 3D
            resadj = dropdims(res, dims=Tuple(ndims(res)-sizediff+1:ndims(res)))
            @test size(resadj) == size(exp)
            @test resadj ≈ exp
        end

        function test_named_graph(g_org, extradims = (); serialize_insizes=false)
            outsize = nout(g_org.inputs[1])
            bs = 4
            indata = reshape(collect(Float32, 1:outsize*bs*prod(extradims)), extradims..., outsize, :)

            gp_org = if serialize_insizes
                graphproto(g_org, name(inputs(g_org)[]) => size(indata); namestrat=ONNXNaiveNASflux.default_namestrat(g_org))
            else
                graphproto(g_org)
            end

            gp_org.name="testmodel"
            validate(modelproto(;graph=gp_org))
            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)

            @test name.(vertices(g_org)) == name.(vertices(g_new))
            @test nout.(vertices(g_org)) == nout.(vertices(g_new))
            @test nin.(vertices(g_org)) == nin.(vertices(g_new))

            expout = g_org(indata)
            resout = g_new(indata)

            test_outputs(resout, expout)

            # For FLux recurrent layers as they accept 2D input but ONNX wants 3D input
            sizediff = length(ONNXNaiveNASflux.shape(g_new.inputs[])) - ndims(indata)
            indata = reshape(indata, size(indata)..., ones(Int, sizediff)...)

            ortout = onnxruntime_infer(g_org, indata)

            test_outputs(ortout, expout, sizediff)

            return g_new
        end

        @testset "Generic function" begin
            l1 = Dense(2, 3, elu)
            l2 = Dense(3, 2)
            function f(x, y)
                x1 = l1(x)
                x2 = l2(x1)
                return x2 .+ y
            end

            gp_sizes = graphproto(f, "x" => (2, missing), "y" => (2,))
            g_sizes = CompGraph(serdeser(gp_sizes))

            x = reshape(collect(Float32, 1:2*4), 2,4)
            y = Float32[5, 6]

            @test name.(vertices(g_sizes)) == ["x", "dense_0", "dense_1", "y", "add_0"]
            @test nout.(vertices(g_sizes)) == [2, 3, 2, 2, 2]
            @test nin.(vertices(g_sizes)) == [[], [2], [3], [], [2, 2]]

            @test g_sizes(x, y) ≈ f(x,y)

            function f(x)
                x1 = l1(x)
                return l2(x1)
            end

            gp_nosizes =  graphproto(f, "x" => missing)
            g_nosizes = CompGraph(serdeser(gp_nosizes))

            @test name.(vertices(g_nosizes)) == ["x", "dense_0", "dense_1"]
            @test g_nosizes(x) ≈ f(x)
        end

        @testset "Linear Dense graph" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 4, relu)
            v2 = dense("dense2", v1, 5, elu)
            v3 = dense("output", v2, 2)

            test_named_graph(CompGraph(v0, v3))
        end

        @testset "Linear Dense graph without names" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense(v0, 4, selu)
            v2 = dense(v1, 5, relu)
            v3 = dense(v2, 2)

            g_org = CompGraph(v0, v3)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)

            @test name.(vertices(g_new)) == ["input_0", "dense_0", "dense_1", "dense_2"]

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "Linear Dense graph non-unique names" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense("vv", v0, 4, selu)
            v2 = dense("vv", v1, 5, relu)
            g_org = CompGraph(v0, v2)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)

            @test name.(vertices(g_new)) == ["input_0", "dense_0", "dense_1",]

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "Linear Conv graph" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv1", v0, 4, selu)
            v2 = convvertex("conv2", v1, 5, elu)
            v3 = convvertex("output", v2, 2)

            test_named_graph(CompGraph(v0, v3), (2,3))
        end

        @testset "Linear Conv graph with global pooling" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv1", v0, 4, relu)
            v2 = convvertex("conv2", v1, 5, elu)
            v3 = fvertex("globmeanpool", v2, x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))))
            v4 = dense("output", v3, 2)

            test_named_graph(CompGraph(v0, v4), (2,3))
        end

        @testset "Linear Conv graph with global pooling without names" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("", v0, 4, relu)
            v2 = invariantvertex(x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))), v1)

            g_org = CompGraph(v0, v2)

            gp_org = graphproto(g_org)
            @test length(size(gp_org.output[])) == 2

            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)
            @test name.(vertices(g_new)) == ["input_0", "conv_0", "globalaveragepool_0"]

            indata = reshape(collect(Float32, 1:3*2*2*2), 2,2,3,2)
            @test size(g_org(indata)) == size(g_new(indata))
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "Linear Batchnorm and Conv graph with global pooling" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv", v0, 4, relu)
            v2 = bnvertex("batchnorm", v1, elu)
            v3 = fvertex("globmeanpool", v2, x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))))
            v4 = dense("output", v3, 2, selu)

            test_named_graph(CompGraph(v0, v4), (4,6))
        end

        @testset "Linear Conv and MaxPool graph with global pooling" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = mpvertex("maxpool", v0)
            v2 = convvertex("conv", v1, 4, relu)
            v3 = fvertex("globmeanpool", v2, x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))))
            v4 = dense("output", v3, 2, selu)

            test_named_graph(CompGraph(v0, v4), (2,3))
        end

        @testset "Dense graph with add" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 4, relu)
            v2 = dense("dense2", v0, 4)
            v3= "add" >> v1 + v2
            v4 = dense("output", v3, 2, selu)

            test_named_graph(CompGraph(v0, v4))
        end

        @testset "Dense graph with add without names" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense(v0, 4, relu)
            v2 = dense(v0, 4)
            v3 = v1 + v2

            g_org = CompGraph(v0, v3)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)

            @test name.(vertices(g_new)) == ["input_0", "dense_0", "dense_1", "add_0"]

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "Dense graph with add and layerfun" begin
            import ONNXNaiveNASflux: create_vertex_default
            v0 = denseinputvertex("input", 3)
            v1 = dense(v0, 4, relu)
            v2 = dense(v0, 4)
            v3 = v1 + v2

            g_org = CompGraph(v0, v3)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            callcnt = 0
            struct CntSpy <: AbstractMutableComp
                f
            end
            function (c::CntSpy)(x...)
                callcnt += 1
                return c.f(x...)
            end
            NaiveNASflux.wrapped(c::CntSpy) = c.f
            NaiveNASflux.layer(c::CntSpy) = layer(c.f)
            NaiveNASflux.layertype(c::CntSpy) = layertype(c.f)

            g_new = CompGraph(gt_new; vfun = (args...) -> create_vertex_default(args...;layerfun=CntSpy))

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            outdata = ones(Float32, nout(v3), size(indata, 2))

            Flux.train!((x,y) -> Flux.mse(g_new(x), y), params(g_new), [(indata, outdata)], Flux.Descent(0.6))
            @test callcnt == nvertices(g_new) - 1
        end

        @testset "Dense graph with cat" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 4, elu)
            v2 = dense("dense2", v0, 4)
            v3 = concat("conc", v1, v2)
            v4 = dense("output", v3, 2, relu)

            test_named_graph(CompGraph(v0, v4))
        end

        @testset "Dense graph with cat and layerfun" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 4, elu)
            v2 = dense("dense2", v0, 4)
            v3 = concat("conc", v1, v2)
            v4 = dense("output", v3, 2, relu)

            g_org = CompGraph(v0, v4)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            callcnt = 0
            struct CntSpy <: AbstractMutableComp
                f
            end
            function (c::CntSpy)(x...)
                callcnt += 1
                return c.f(x...)
            end
            NaiveNASflux.wrapped(c::CntSpy) = c.f
            NaiveNASflux.layer(c::CntSpy) = layer(c.f)
            NaiveNASflux.layertype(c::CntSpy) = layertype(c.f)

            g_new = CompGraph(gt_new; vfun = (args...) -> create_vertex_default(args...;layerfun=CntSpy))

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            outdata = ones(Float32, nout(v4), size(indata, 2))

            Flux.train!((x,y) -> Flux.mse(g_new(x), y), params(g_new), [(indata, outdata)], Flux.Descent(0.6))
            @test callcnt == nvertices(g_new) - 1
        end

        @testset "Dense graph with constant scalar op $op" for op in (+, *)
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 4, elu)
            v2 = fvertex("scale", v1, x -> op.(0.5f0, x))
            v3 = dense("output", v2, 2, relu)

            test_named_graph(CompGraph(v0, v3))
        end

        @testset "Dense graph with constant elemwise array op $op" for op in (+, *)
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense1", v0, 2, elu)
            v2 = fvertex("scale", v1, x -> op.(Float32[0.5, 0.1], x))
            v3 = dense("output", v2, 3, relu)

            test_named_graph(CompGraph(v0, v3))
        end

        @testset "Dense graph with free constant" begin
            import ONNXNaiveNASflux: sourcevertex_with_outputs
            v0 = denseinputvertex("input", 3)
            v1 = dense("dense", v0, 2, elu)
            # Constant must be array for generic check to pass as ONNX opset only supports constant tensors
            # Name is not preserved when serializing SourceVertex. We just sneakily set it to what the autoselected name will be
            v2 =  sourcevertex_with_outputs([213], "constant")

            test_named_graph(CompGraph([v0], [v1,v2]))
        end

        @testset "Conv and batchnorm graph with cat" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv", v0, 2, elu)
            v2 = bnvertex("batchnorm", v0)
            v3 = concat("conc", v1, v2)
            v4 = fvertex("globmeanpool", v3, x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))))
            v5 = dense("output", v4, 2, relu)

            test_named_graph(CompGraph(v0, v5), (2,3))
        end

        @testset "Dense graph with cat without names" begin
            v0 = denseinputvertex("input", 3)
            v1 = dense(v0, 4, relu)
            v2 = dense(v0, 4)
            v3 = concat("conc", v1, v2)

            g_org = CompGraph(v0, v3)

            gp_org = graphproto(g_org)
            gt_new = serdeser(gp_org)

            g_new = CompGraph(gt_new)

            @test name.(vertices(g_new)) == ["input_0", "dense_0", "dense_1", "concat_0"]

            indata = reshape(collect(Float32, 1:3*4), nout(v0), :)
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "RNN to LSTM" begin
            v0 = rnninputvertex("input", 3)
            v1 = fluxvertex("rnn", RNN(nout(v0), 4), v0)
            v2 = fluxvertex("lstm", LSTM(nout(v1), 5), v1)

            test_named_graph(CompGraph(v0, v2))
        end

        @testset "Recurrent to Dense" begin
            v0 = rnninputvertex("input", 3)
            v1 = fluxvertex("rnn", RNN(nout(v0), 4), v0)
            v2 = fluxvertex("lstm", LSTM(nout(v1), 5), v1)
            v3 = dense("dense", v2, 6, elu)

            g_org = CompGraph(v0, v3)
            g_new = CompGraph(serdeser(graphproto(g_org)))

            @test name.(vertices(g_new)) == name.(vertices(g_org))

            indata = reshape(collect(Float32, 1:3*5*7), 3,5,7)

            expout = g_org.(Flux.unstack(indata, 3))
            resout = g_new.(Flux.unstack(indata, 3))

            @test size.(expout) == size.(resout)
            @test expout ≈ resout

            ortout, = onnxruntime_infer(g_org, indata)
            expout_s = hcat(expout...)

            @test size(expout_s) == size(ortout)
            @test expout_s ≈ ortout
        end

        @testset "Graph two inputs two outputs" begin
            vins = denseinputvertex.(["in1", "in2"], 3)
            v1 = "add" >> vins[1] + vins[2]
            v2 = concat("conc", vins[1], vins[2])

            g_org = CompGraph(vins, [v1, v2])
            g_new = CompGraph(serdeser(graphproto(g_org)))

            @test name.(vertices(g_org)) == name.(vertices(g_new))

            indata1 = reshape(collect(Float32, 1:3*4), nout(vins[1]), :)
            indata2 = indata1 .* -0.5
            @test g_org(indata1, indata2) == g_new(indata1, indata2)
        end

        @testset "Conv to reshape to Dense" begin
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv", v0, 4, relu)
            v2 = absorbvertex("reshape", ONNXNaiveNASflux.Reshape((12, Colon())), v1; traitdecoration=SizePseudoTransparent)
            v3 = dense("dense", v2, 3, elu)

            test_named_graph(CompGraph(v0, v3), (2,3))
        end

        @testset "Conv to Flatten to Dense" begin
            using ONNXNaiveNASflux: MeasureNout
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv", v0, 4, relu)
            v2 = absorbvertex("flatten",MeasureNout(ONNXNaiveNASflux.Flatten(3); outsize=24), v1; traitdecoration=SizePseudoTransparent)
            v3 = dense("dense", v2, 3, elu)

            test_named_graph(CompGraph(v0, v3), (2,3); serialize_insizes=true)
        end

        @testset "Conv to flatten to Dense" begin
            using ONNXNaiveNASflux: MeasureNout
            v0 = conv2dinputvertex("input", 3)
            v1 = convvertex("conv", v0, 4, relu)
            v2 = absorbvertex("flatten", MeasureNout(Flux.flatten; actdim=1, outsize=24), v1; traitdecoration=SizePseudoTransparent)
            v3 = dense("dense", v2, 3, elu)

            test_named_graph(CompGraph(v0, v3), (2,3); serialize_insizes=true)
        end

        @testset "Infer shapes" begin

            function remodel(g, args...=missing)
                pb = PipeBuffer()
                save(pb, g, args...)
                load(pb)
            end

            @testset "Batchnorm -> Conv graph" begin
                v0 = conv2dinputvertex("input", 3)
                v1 = bnvertex("v1", v0)
                v2 = convvertex("v2", v1, 2)

                g = remodel(CompGraph(v0, v2))
                @test layertype(g.inputs[1]) == layertype(v0)
            end

            @testset "Two conv graph" begin
                v0 = conv2dinputvertex("input", 3)
                v1a = convvertex("v1a", v0, 4)
                v1b = convvertex("v1b", v0, 2)
                v2 = concat("v2", v1a, v1b)

                g_org = g = remodel(CompGraph(v0, v2))
                @test layertype(g.inputs[1]) == layertype(v0)
            end

            @testset "Concat path" begin
                v0 = denseinputvertex("input", 3)
                v1 = dense("v1", v0, 4, elu)
                v2 = concat(v1, v0)

                g = remodel(CompGraph(v0, v2))
                @test layertype(g.inputs[1]) == layertype(v0)
            end

            @testset "Shortcut to globpool -> dense" begin
                v0 = conv2dinputvertex("input", 3)
                v1 = convvertex("v1", v0, 2)
                v2 = concat("v2", v1, v0)
                v3 = fvertex("v3", v2, x -> ONNXNaiveNASflux.globalmeanpool(x, y -> dropdims(y, dims=(1,2))))
                v4 = dense("v4", v3, 4)

                g = remodel(CompGraph(v0, v4))
                @test layertype(g.inputs[1]) == layertype(v0)
            end
        end
    end

    @testset "Models" begin
        import ONNXNaiveNASflux: modelproto, sizes

        @testset "Generic function infer" begin
            _f(x, y) = x .+ y
            f(x, y) = _f(x,y)
            f(x::Matrix{Int}, y::Matrix{Int}) = _f(x, y)

            mp = modelproto(f)
            mt = serdeser(mp)
            ss = sizes(mt)

            @test length(ss["data_0"]) == 2
            @test length(ss["data_1"]) == 2

            g = @test_logs (:warn, r"No valid input sizes provided.") load(mt)
            g([1,2], [3,4]) == f([1,2], [3,4])
        end

        @testset "CompGraph infer" begin
            vis = denseinputvertex.("in", [2, 2])
            g_org = CompGraph(vis, +(vis...))

            mp = modelproto(g_org)
            mt = serdeser(mp)
            ss = sizes(mt)

            @test length(ss["in_0"]) == 2
            @test length(ss["in_1"]) == 2

            g_new = load(mt)
            g_org([1,2], [3,4]) == g_new([1,2], [3,4])
        end

        @testset "$(cfun(l)) infer" for (l, expshape) in (
            (Dense(2,3), (2,0)),
            (Conv((1, 1), 2=>3), (0,0,2,0)),
            (RNN(2,3), (2,0,0)),
            (LSTM(2,3), (2,0,0)),
            (SkipConnection(Dense(2,2), +), (2,0)),
            ), cfun in (Chain, l -> Chain(BatchNorm(nin(l)[]), l))

            c_org = cfun(l)

            mt = modelproto(c_org) |> serdeser
            ss = sizes(mt)

            @test ss["data_0"] == expshape
        end
    end

    @testset "Save to file" begin
        using ONNXNaiveNASflux.NaiveNASflux
        function tryfile(filename, args...; kwargs...)
            try
                save(filename, args...; kwargs...)
                return load(filename)
            finally
                rm(filename;force=true)
            end
        end

        @testset "Generic function no pars" begin
            f = (x,y) -> x + y
            g = tryfile("generic_function_no_pars.onnx", f)
            @test name.(vertices(g)) == ["data_0", "data_1", "add_0"]
            @test g(1,3) == f(1, 3)
        end

        @testset "Generic function sizes" begin
            f = (x,y) -> x + y
            g = tryfile("generic_function_sizes.onnx", f, (2,missing), (2,missing))
            nout.(vertices(g)) == [2, 2]
            @test g([1,3], [2, 4]) == f([1,3], [2, 4])
        end

        @testset "Simple graph no pars" begin
            v0 = denseinputvertex("in", 3)
            v1 = fluxvertex("dense1", Dense(3, 2, relu), v0)
            v2 = fluxvertex("dense2", Dense(2, 3), v1)
            g_org = CompGraph(v0, v2)

            g_new = tryfile("simple_graph_no_pars.onnx", g_org)

            @test name.(vertices(g_org)) == name.(vertices(g_new))
            @test nout.(vertices(g_org)) == nout.(vertices(g_new))
            @test nin.(vertices(g_org)) == nin.(vertices(g_new))

            indata = reshape(collect(Float32, 1:3*2), 3,2)
            @test g_org(indata) ≈ g_new(indata)
        end

        @testset "Simple graph namestrat" begin
            v0 = conv2dinputvertex("in", 3)
            v1 = fluxvertex("conv", Conv((1,2), 3 => 4, relu), v0)
            v2 = fluxvertex("bn", BatchNorm(4, elu), v1)
            g_org = CompGraph(v0, v2)

            ng = ONNXNaiveNASflux.name_runningnr()
            ns(::NaiveNASlib.MutationVertex) = n -> ng
            ns(n) = ng(n)

            g_new = tryfile("simple_graph_namestrat.onnx", g_org; namestrat=ns)

            @test name.(vertices(g_new)) == ["in_0", "conv_0", "batchnorm_0"]
            @test nout.(vertices(g_org)) == nout.(vertices(g_new))
            @test nin.(vertices(g_org)) == nin.(vertices(g_new))

            indata = reshape(collect(Float32, 1:1*2*3*4), 1,2,3,4)
            @test g_org(indata) ≈ g_new(indata)
        end
    end

    @testset "Allowed input shapes" begin
        function remodel(m, args...; assertwarn=true)
            pb = PipeBuffer()
            save(pb, m, args...)
            if assertwarn
                return @test_logs (:warn, r"No valid input sizes") load(pb)
            end
            return load(pb)
        end

        @testset "Allowed input shapes op: $(tc[1])" for tc in (
            (Dense(2, 3), ((2, 3), (2, missing), missing, (missing, missing)), (2, 2)),
            (Conv((1,), 2=>3), ((1,2,1), (1, missing, missing), missing, ntuple(i -> missing, 3)), (1,2,1)),
            (Conv((1,1), 2=>3), ((1,1,2,1), (1, missing, missing, missing), missing, ntuple(i -> missing, 4)), (1,1,2,1)),
            (RNN(2,3), ((2,3,1), missing, ntuple(i -> missing, 3), ntuple(i -> missing ,4)), (2, 3))
        )
            op, testsizes, validsize = tc
            inpt = ones(Float32, validsize)

            g1 = remodel(op; assertwarn=false)
            @test g1(inpt) == op(inpt)
            @testset "Inputshape $s" for s in testsizes
                assertwarn = s isa Tuple && length(s) != length(shape(layertype(op), 1))
                g = remodel(op, s; assertwarn)
                Flux.reset!(op) # For RNNs or else the test below will fail
                @test g(inpt) == op(inpt)
            end
        end

        @testset "Allowed input shapes op: +" begin
            in1, in2 = ones(3), ones(3)
            op = (x,y) -> x .+ y
            g1 = remodel(op; assertwarn=false)
            @test g1(in1, in2) == in1 .+ in2
            @testset "Inputshape $s1, $s2" for (s1, s2, assertwarn) in (
                ((3,), (3,), false),
                ((missing,), (missing,), true),
                (missing, missing, false),
                )
                g = remodel(op, s1, s2; assertwarn)
                @test g(in1, in2) == in1 .+ in2
            end
        end

        @testset "Disallowed input shapes op: $(tc[1])" for tc in (
            (Dense(2,3), ((2,), (3, 1), (missing, ), (1,2,3,4))),
            (Conv((1,), 2=>3), ((2,), (missing, ), (1,1,1), (2,3,4,5,6), (2,3,4,5))),
            (Conv((1,1), 2=>3), ((2,), (missing, ), (1,1,1,1), (2,3,4,5,6), (2,3,4))),
            (MaxPool((2,2)), ((2,), (missing, ), (1,1,1), (2,3,4,5,6))),
            (RNN(2,3), ((2, ), (missing,), (2,1), (1,2,3), (2,3,4,5,6)))
        )
            op, testsizes = tc
            @testset "Inputshape $s" for s in testsizes
                @test_throws DimensionMismatch remodel(op, s)
            end
        end
    end
end


text:
const creamery = RGB(214/255, 210/255, 196/255)
const stegeman = RGB(158/255, 162/255, 162/255)

H_radius = Node{Float32}(0.3)

function get_atom_radius(S::String)
    if S == "H"
        return to_value(H_radius)
    else
        return 0.5
    end
end

function get_atom_color(S::String)
    if haskey(atom_color, S)
        return atom_color[S]
    else
        throw(ArgumentError("Color not defined for atom $S"))
    end
end

const van_der_waals = Dict(
    "H" => 32,
    "He" => 46,
    "Li" => 133,
    "Be" => 102,
    "B" => 85,
    "C" => 75,
    "N" => 71,
    "O" => 63,
    "F" => 64,
    "Ne" => 67,
    "Na" => 155,
    "Mg" => 139,
    "Al" => 126,
    "Si" => 116,
    "P" => 111,
    "S" => 103,
    "Cl" => 99,
    "Ar" => 96,
    "K" => 196,
    "Ca" => 171,
    "Sc" => 148,
    "Ti" => 136,
    "V" => 134,
    "Cr" => 122,
    "Mn" => 120,
    "Fe" => 120,
    "Co" => 120,
    "Ni" => 120,
    "Cu" => 120,
    "Zn" => 120,
    "Ga" => 120,
    "Ge" => 120,
    "As" => 120,
    "Se" => 120,
    "Br" => 120,
    "Kr" => 120,
    "Rb" => 120,
    "Sr" => 120,
    "Y" => 120,
    "Zr" => 120,
    "Nb" => 120,
    "Mo" => 120,
    "Tc" => 120,
    "Ru" => 120,
    "Rh" => 120,
    "Pd" => 120,
    "Ag" => 120,
    "Cd" => 120,
    "In" => 120,
    "Sn" => 120,
    "Sb" => 120,
    "Te" => 120,
    "I" => 120,
    "Xe" => 120,
    "Cs" => 120,
    "Ba" => 120,
    "La" => 120,
    "Ce" => 120,
    "Pr" => 120,
    "Nd" => 120,
    "Pm" => 120,
    "Sm" => 120,
    "Eu" => 120,
    "Gd" => 120,
    "Tb" => 120,
    "Dy" => 120,
    "Ho" => 120,
    "Er" => 120,
    "Tm" => 120,
    "Yb" => 120,
    "Lu" => 120,
    "Hf" => 120,
    "Ta" => 120,
    "W" => 120,
    "Re" => 120,
    "Os" => 120,
    "Ir" => 120,
    "Pt" => 120,
    "Au" => 120,
    "Hg" => 120,
    "Tl" => 120,
    "Pb" => 120,
    "Bi" => 120,
    "Po" => 120,
    "At" => 120,
    "Rn" => 120,
    "Fr" => 120,
    "Ra" => 120,
    "Ac" => 120,
    "Th" => 120,
    "Pa" => 120,
    "U" => 120,
    "Np" => 120,
    "Pu" => 120,
    "Am" => 120,
    "Cm" => 120,
    "Bk" => 120,
    "Cf" => 120,
    "Es" => 120,
    "Fm" => 120,
    "Md" => 120,
    "No" => 120,
    "Lr" => 120,
    "Rf" => 120,
    "Db" => 120,
    "Sg" => 120,
    "Bh" => 120,
    "Hs" => 120,
    "Mt" => 120,
    "Ds" => 120,
    "Rg" => 120,
    "Cn" => 120,
    "Nh" => 120,
    "Fl" => 120,
    "Mc" => 120,
    "Lv" => 120,
    "Ts" => 120,
    "Og" => 120
)

const atom_color = Dict(
    "H" => color("whitesmoke"),
    "He" => 46,
    "Li" => 133,
    "Be" => 102,
    "B" => 85,
    "C" => color("gray19"),
    "N" => color("royalblue1"),
    "O" => RGB(0.7294117647058823,0.047058823529411764,0.1843137254901961),
    "F" => 64,
    "Ne" => 67,
    "Na" => 155,
    "Mg" => 139,
    "Al" => 126,
    "Si" => 116,
    "P" => 111,
    "S" => color("darkorange"),
    "Cl" => 99,
    "Ar" => 96,
    "K" => 196,
    "Ca" => 171
)

ext: jl

max_stars_repo_name: gustavojra/Alchemy

max_issues_repo_name: gustavojra/Alchemy

max_forks_repo_name: gustavojra/Alchemy


text:
struct RectilinearMeshPartition
    grid::RectilinearGrid{X,Y,Z,T}

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testset "6_serialization" begin
  parms_stream = UInt8[]
  data_stream1 = UInt8[]
  data_stream2 = UInt8[]
  data_stream3 = UInt8[]
  data_stream4 = UInt8[]
  sk_stream = UInt8[]

  @testset "server (part 1)" begin
    enc_parms = EncryptionParameters(SchemeType.ckks)
    poly_modulus_degree = 8192
    @testset "polynomial modulus degree" begin
      @test_nowarn set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
    end

    @testset "coefficient modulus" begin
      @test_nowarn set_coeff_modulus!(enc_parms, coeff_modulus_create(poly_modulus_degree, [50, 20, 50]))
    end

    @testset "save! EncryptionParameters" begin
      @test save_size(enc_parms) == 192
      resize!(parms_stream, save_size(enc_parms))
      @test save!(parms_stream, enc_parms) == 75
      out_bytes = 75
      resize!(parms_stream, out_bytes)
    end

    @testset "save_size comparison" begin
      @test save_size(ComprModeType.none, enc_parms) == 129
      @test save_size(ComprModeType.zlib, enc_parms) == 146
    end

    @testset "save! and load! EncryptionParameters" begin
      byte_buffer = Vector{UInt8}(undef, save_size(enc_parms))
      @test save!(byte_buffer, length(byte_buffer), enc_parms) == 75

      enc_parms2 = EncryptionParameters()
      @test load!(enc_parms2, byte_buffer, length(byte_buffer)) == 75
      @test enc_parms == enc_parms2
    end
  end

  @testset "client (part 1)" begin
    enc_parms = EncryptionParameters()
    @testset "load! EncryptionParameters" begin
      @test load!(enc_parms, parms_stream) == 75
    end

    context = SEALContext(enc_parms)
    keygen = KeyGenerator(context)
    pk = PublicKey()
    create_public_key!(pk, keygen)
    sk = secret_key(keygen)

    @testset "save! SecretKey" begin
      @test save_size(sk) == 197464
      resize!(sk_stream, save_size(sk))
      @test isapprox(save!(sk_stream, sk), 145823, rtol=0.01)
      out_bytes = save!(sk_stream, sk)
      resize!(sk_stream, out_bytes)
    end

    rlk = create_relin_keys(keygen)
    @testset "save! create_relin_keys" begin
      @test save_size(rlk) == 395189
      resize!(data_stream1, save_size(rlk))
      @test isapprox(save!(data_stream1, rlk), 291635, rtol=0.01)
      size_rlk = save!(data_stream1, rlk)
      resize!(data_stream1, size_rlk)
    end

    rlk_big = RelinKeys()
    create_relin_keys!(rlk_big, keygen)
    @testset "save! create_relin_keys" begin
      @test save_size(rlk_big) == 789779
      resize!(data_stream2, save_size(rlk_big))
      @test isapprox(save!(data_stream2, rlk_big), 583244, rtol=0.01)
      size_rlk_big = save!(data_stream2, rlk_big)
      resize!(data_stream2, size_rlk_big)
    end

    initial_scale = 2.0^20
    encoder = CKKSEncoder(context)
    plain1 = Plaintext()
    plain2 = Plaintext()

    @testset "encode!" begin
      @test_nowarn encode!(plain1, 2.3, initial_scale, encoder)
      @test_nowarn encode!(plain2, 4.5, initial_scale, encoder)
    end

    encryptor = Encryptor(context, pk)
    encrypted1 = Ciphertext()
    @testset "encrypt!" begin
      @test_nowarn encrypt!(encrypted1, plain1, encryptor)
    end

    @testset "set_secret_key!" begin
      @test_nowarn set_secret_key!(encryptor, sk)
    end

    @testset "encrypt_symmetric, encrypt_symmetric!" begin
      @test encrypt_symmetric(plain2, encryptor) isa Ciphertext
      c = Ciphertext()
      @test encrypt_symmetric!(c, plain2, encryptor) == c
    end
    sym_encrypted2 = encrypt_symmetric(plain2, encryptor)

    @testset "save! Ciphertext" begin
      @test save_size(encrypted1) == 263273
      resize!(data_stream2, save_size(encrypted1))
      @test isapprox(save!(data_stream2, encrypted1), 173531, rtol=0.01)
      size_encrypted1 = save!(data_stream2, encrypted1)
      resize!(data_stream2, size_encrypted1)

      @test save_size(sym_encrypted2) == 131770
      resize!(data_stream3, save_size(sym_encrypted2))
      @test isapprox(save!(data_stream3, sym_encrypted2), 86966, rtol=0.01)
      size_sym_encrypted2 = save!(data_stream3, sym_encrypted2)
      resize!(data_stream3, size_sym_encrypted2)
    end
  end

  @testset "server (part 2)" begin
    enc_parms = EncryptionParameters()
    @testset "load! EncryptionParameters" begin
      @test load!(enc_parms, parms_stream) == 75
    end

    context = SEALContext(enc_parms)
    evaluator = Evaluator(context)
    rlk = RelinKeys()
    encrypted1 = Ciphertext()
    encrypted2 = Ciphertext()

    @testset "load! RelinKeys" begin
      @test isapprox(load!(rlk, context, data_stream1), 291635, rtol=0.01)
    end

    @testset "load! Ciphertext" begin
      @test isapprox(load!(encrypted1, context, data_stream2), 173702, rtol=0.01)
      @test isapprox(load!(encrypted2, context, data_stream3), 86966, rtol=0.01)
    end

    encrypted_prod = Ciphertext()
    @testset "multiply, relinearize, rescale" begin
      @test multiply!(encrypted_prod, encrypted1, encrypted2, evaluator) == encrypted_prod
      @test relinearize_inplace!(encrypted_prod, rlk, evaluator) == encrypted_prod
      @test rescale_to_next_inplace!(encrypted_prod, evaluator) == encrypted_prod
    end

    @testset "save! Ciphertext" begin
      @test save_size(encrypted_prod) == 131689
      resize!(data_stream4, save_size(encrypted_prod))
      @test isapprox(save!(data_stream4, encrypted_prod), 117909, rtol=0.01)
      size_encrypted_prod = save!(data_stream4, encrypted_prod)
      resize!(data_stream4, size_encrypted_prod)
    end
  end

  @testset "client (part 2)" begin
    enc_parms = EncryptionParameters()
    load!(enc_parms, parms_stream)
    context = SEALContext(enc_parms)

    sk = SecretKey()
    @testset "load! SecretKey" begin
      @test isapprox(load!(sk, context, sk_stream), 145823, rtol=0.01)
    end

    decryptor = Decryptor(context, sk)
    encoder = CKKSEncoder(context)
    encrypted_result = Ciphertext()
    @testset "load! Ciphertext" begin
      @test_nowarn load!(encrypted_result, context, data_stream4)
    end

    plain_result = Plaintext()
    @testset "decrypt!" begin
      @test_nowarn decrypt!(plain_result, encrypted_result, decryptor)
    end

    slot_count_ = slot_count(encoder)
    result = Vector{Float64}(undef, slot_count_)
    @testset "decode! and check result" begin
      @test_nowarn decode!(result, plain_result, encoder)
      @test isapprox(result[1], 10.35, rtol=0.01)
      @test isapprox(result[2], 10.35, rtol=0.01)
      @test isapprox(result[3], 10.35, rtol=0.01)
      @test isapprox(result[end-2], 10.35, rtol=0.01)
      @test isapprox(result[end-1], 10.35, rtol=0.01)
      @test isapprox(result[end-0], 10.35, rtol=0.01)
    end
  end

  pt = Plaintext("1x^2 + 3")
  stream = Vector{UInt8}(undef, save_size(pt))
  @testset "save! Plaintext" begin
    @test save!(stream, pt) == 66
    data_size = 66
    resize!(stream, data_size)
  end

  header = SEALHeader()
  @testset "load_header!" begin
    @test load_header!(header, stream) == header
    @test header.size == 66
  end
end


text:
"""
    SerialExecution

Functor object for using serial execution in a lshaped algorithm. Create by supplying a [`Serial`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SerialExecution{H <: AbstractFeasibilityHandler,
                       T <: AbstractFloat,
                       A <: AbstractVector} <: AbstractLShapedExecution
    subproblems::Vector{SubProblem{H,T}}
    decisions::Decisions
    subobjectives::A
    model_objectives::A

    function SerialExecution(structure::VerticalStructure{2, 1, <:Tuple{ScenarioProblems}},
                             ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                     T <: AbstractFloat,
                                                                     A <: AbstractVector}
        H = HandlerType(F)
        return new{H,T,A}(Vector{SubProblem{H,T}}(), structure.decisions[2], A(), A())
    end
end

function num_thetas(lshaped::AbstractLShaped, ::SerialExecution)
    return num_thetas(num_subproblems(lshaped), lshaped.aggregation)
end

function initialize_subproblems!(execution::SerialExecution{H,T},
                                 scenarioproblems::ScenarioProblems) where {H <: AbstractFeasibilityHandler,
                                                                            T <: AbstractFloat}
    for i in 1:num_subproblems(scenarioproblems)
        push!(execution.subproblems, SubProblem(
            subproblem(scenarioproblems, i),
            i,
            T(probability(scenario(scenarioproblems, i))),
            execution.decisions.knowns,
            H))
    end
    return nothing
end

function finish_initilization!(lshaped::AbstractLShaped, execution::SerialExecution)
    append!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
    append!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
    return nothing
end

function restore_subproblems!(lshaped::AbstractLShaped, execution::SerialExecution)
    for subproblem in execution.subproblems
        restore_subproblem!(subproblem)
    end
    return nothing
end

function resolve_subproblems!(lshaped::AbstractLShaped, execution::SerialExecution{H,T}) where {H <: AbstractFeasibilityHandler, T <: AbstractFloat}
    # Update subproblems
    update_known_decisions!(execution.decisions, lshaped.x)
    change = KnownValuesChange()
    # Assume no cuts are added
    added = false
    # Update and solve subproblems
    for subproblem in execution.subproblems
        update_subproblem!(subproblem, change)
        cut::SparseHyperPlane{T} = subproblem(lshaped.x)
        added |= aggregate_cut!(lshaped, lshaped.aggregation, cut)
    end
    added |= flush!(lshaped, lshaped.aggregation)
    # Return current objective value and cut_added flag
    return current_objective_value(lshaped), added
end

# API
# ------------------------------------------------------------
function (execution::Serial)(structure::VerticalStructure{2, 1, <:Tuple{ScenarioProblems}},
                             ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                     T <: AbstractFloat,
                                                                     A <: AbstractVector}
    return SerialExecution(structure, F, T, A)
end

function str(::Serial)
    return ""
end


text:
export
    Packed_Factor

mutable struct Packed_Factor
    measurement::Vector{Dict{String, Any}}
    distType::String
end


text:
function assert_option(x)
    is_option(x) || error("$(typeof(x)) is not an option type")
end

@option struct ConvertOption
    include_defaults::Bool=false
    exclude_nothing::Bool=false
end

"""
    to_dict(option; include_defaults=false, exclude_nothing=false) -> OrderedDict

Convert an object to an `OrderedDict`.

# Kwargs

- `include_defaults`: include the default value.
- `exclude_nothing`: exclude fields that have value `nothing`,
    this supersedes `include_defaults` when they are both `true`.

!!! tips
    `to_dict` does not export fields that are of the same values as the defaults.
    In most cases, this should be the default behaviour, and users should not use
    `include_defaults`, however,  this can be overridden by changing `include_defaults`
    to `true`.
"""
function to_dict(x; include_defaults::Bool=false, exclude_nothing::Bool=false)
    assert_option(x)
    return to_dict(typeof(x), x, ConvertOption(include_defaults, exclude_nothing))
end

"""
    to_dict(::Type{T}, x, option::ConvertOption) where T

Convert `x` when `x` is inside an option type `T`. `option`
is a set of options to determine the conversion behaviour. this can
be overloaded to change the behaviour of `to_dict(x; kw...)`.

    to_dict(::Type{T}, x) where T

One can also use the 2-arg version when `x` is not or does not
contain an option type for convenience.

# Example

The following is a builtin overload to handle list of options.

```julia
function Configurations.to_dict(::Type{T}, x::Vector, option::ConvertOption) where T
    if is_option(eltype(x))
        return map(p->to_dict(T, p, include_defaults), x)
    else
        return x
    end
end
```

The following overloads the 2-arg `to_dict` to convert all `VersionNumber` to
a `String` for all kinds of option types.

```julia
Configurations.to_dict(::Type, x::VersionNumber) = string(x)
```
"""
function to_dict(::Type{T}, x, option::ConvertOption) where T
    if is_option(x)
        return _option_to_dict(x, option)
    else
        return to_dict(T, x) # fall through 2-arg version
    end
end

to_dict(::Type, x) = x

# handle list of options as builtin
function to_dict(::Type{T}, x::Vector, option::ConvertOption) where T
    if is_option(eltype(x))
        return map(p->to_dict(T, p, option), x)
    else
        return x
    end
end

function _option_to_dict(x, option::ConvertOption)
    assert_option(x)

    d = OrderedDict{String, Any}()
    T = typeof(x)
    for name in fieldnames(T)
        type = fieldtype(T, name)
        value = getfield(x, name)
        if option.exclude_nothing && value === nothing
            continue
        end

        if option.include_defaults || value != field_default(T, name)
            field_dict = to_dict(T, value, option)

            # always add an alias if it's a Union
            # of multiple option types
            if is_option(value) && type isa Union
                d[string(name)] = OrderedDict{String, Any}(
                    type_alias(typeof(value)) => field_dict,
                )
            else
                d[string(name)] = field_dict
            end
        end
    end
    return d
end

"""
    to_toml([f::Function], io::IO, option; sorted=false, by=identity, kw...)

Convert an instance `option` of option type to TOML and write it to `IO`. See [`to_dict`](@ref)
for other valid keyword options. See also `TOML.print` in the stdlib for the explaination of
`sorted`, `by` and `f`.
"""
function to_toml(f, io::IO, x; sorted::Bool=false, by=identity, kw...)
    is_option(x) || error("argument is not an option type")
    d = to_dict(x; kw...)
    return TOML.print(f, io, d; sorted=sorted, by=by)
end

"""
    to_toml([f::Function], filename::String, option; sorted=false, by=identity, kw...)

Convert an instance `option` of option type to TOML and write it to `filename`. See also `TOML.print`.
"""
function to_toml(f, filename::String, x; sorted::Bool=false, by=identity, kw...)
    open(filename, "w+") do io
        to_toml(f, io, x; sorted=sorted, by=by, kw...)
    end
end

function to_toml(io::IO, x; sorted::Bool=false, by=identity, kw...)
    return to_toml(identity, io, x; sorted=sorted, by=by, kw...)
end

function to_toml(filename::String, x; sorted::Bool=false, by=identity, kw...)
    return to_toml(identity, filename, x; sorted=sorted, by=by, kw...)
end

"""
    to_toml(x; sorted=false, by=identity, kw...)

Convert an instance `x` of option type to TOML and write it to `String`. See also `TOML.print`.

`to_toml` does not export fields that are of the same values as the defaults. This can be
overridden by changing `include_defaults` to `true`.
"""
function to_toml(x; sorted::Bool=false, by=identity, kw...)
    return sprint(x) do io, x
        to_toml(io, x; sorted=sorted, by=by, kw...)
    end
end


text:
using ClimateMachine.VariableTemplates: Vars, Grad, @vars
using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws:
    BalanceLaw,
    vars_state,
    state_to_entropy_variables!,
    entropy_variables_to_state!,
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    state_to_entropy,
    boundary_conditions,
    boundary_state!,
    wavespeed,
    flux_first_order!,
    source!

using StaticArrays
using LinearAlgebra: dot, I
using ClimateMachine.DGMethods.NumericalFluxes: NumericalFluxFirstOrder
import ClimateMachine.DGMethods.NumericalFluxes:
    EntropyConservative,
    numerical_volume_conservative_flux_first_order!,
    numerical_volume_fluctuation_flux_first_order!,
    ave,
    logave,
    numerical_flux_first_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_second_order!

using ClimateMachine.Orientations:
    Orientation, FlatOrientation, SphericalOrientation
using ClimateMachine.Atmos: NoReferenceState
using ClimateMachine.Grids

using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet: grav, R_d, cp_d, cv_d, planet_radius, MSLP, Omega

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

const total_energy = false
const fluctuation_gravity = false

@inline gamma(ps::EarthParameterSet) = cp_d(ps) / cv_d(ps)

abstract type AbstractDryAtmosProblem end

struct DryAtmosModel{D, O, P, RS, S, DS} <: BalanceLaw
    orientation::O
    problem::P
    ref_state::RS
    sources::S
    drag_source::DS
end
function DryAtmosModel{D}(
    orientation,
    problem::AbstractDryAtmosProblem;
    ref_state = NoReferenceState(),
    sources = (),
    drag_source = NoDrag(),
) where {D}
    O = typeof(orientation)
    P = typeof(problem)
    RS = typeof(ref_state)
    S = typeof(sources)
    DS = typeof(drag_source)
    DryAtmosModel{D, O, P, RS, S, DS}(
        orientation,
        problem,
        ref_state,
        sources,
        drag_source,
    )
end

boundary_conditions(::DryAtmosModel) = (1, 2)
# XXX: Hack for Impenetrable.
#      This is NOT entropy stable / conservative!!!!
function boundary_state!(
    ::NumericalFluxFirstOrder,
    bctype,
    ::DryAtmosModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    _...,
)
    state⁺.ρ = state⁻.ρ
    state⁺.ρu -= 2 * dot(state⁻.ρu, n) .* SVector(n)
    state⁺.ρe = state⁻.ρe
    aux⁺.Φ = aux⁻.Φ
end

function init_state_prognostic!(m::DryAtmosModel, args...)
    init_state_prognostic!(m, m.problem, args...)
end

function nodal_init_state_auxiliary!(
    m::DryAtmosModel,
    state_auxiliary,
    tmp,
    geom,
)
    init_state_auxiliary!(m, m.orientation, state_auxiliary, geom)
    init_state_auxiliary!(m, m.ref_state, state_auxiliary, geom)
    init_state_auxiliary!(m, m.problem, state_auxiliary, geom)
end

function altitude(::DryAtmosModel{dim}, ::FlatOrientation, geom) where {dim}
    @inbounds geom.coord[dim]
end

function altitude(::DryAtmosModel, ::SphericalOrientation, geom)
    FT = eltype(geom)
    _planet_radius::FT = planet_radius(param_set)
    norm(geom.coord) - _planet_radius
end

"""
    init_state_auxiliary!(
        m::DryAtmosModel,
        aux::Vars,
        geom::LocalGeometry
        )

Initialize geopotential for the `DryAtmosModel`.
"""
function init_state_auxiliary!(
    ::DryAtmosModel{dim},
    ::FlatOrientation,
    state_auxiliary,
    geom,
) where {dim}
    FT = eltype(state_auxiliary)
    _grav = FT(grav(param_set))
    @inbounds r = geom.coord[dim]
    state_auxiliary.Φ = _grav * r
    state_auxiliary.∇Φ =
        dim == 2 ? SVector{3, FT}(0, _grav, 0) : SVector{3, FT}(0, 0, _grav)
end
function init_state_auxiliary!(
    ::DryAtmosModel,
    ::SphericalOrientation,
    state_auxiliary,
    geom,
)
    FT = eltype(state_auxiliary)
    _grav = FT(grav(param_set))
    r = norm(geom.coord)
    state_auxiliary.Φ = _grav * r
    state_auxiliary.∇Φ = _grav * geom.coord / r
end

function init_state_auxiliary!(
    ::DryAtmosModel,
    ::NoReferenceState,
    state_auxiliary,
    geom,
) end

function init_state_auxiliary!(
    ::DryAtmosModel,
    ::AbstractDryAtmosProblem,
    state_auxiliary,
    geom,
) end

struct DryReferenceState{TP}
    temperature_profile::TP
end
vars_state(::DryAtmosModel, ::DryReferenceState, ::Auxiliary, FT) =
    @vars(T::FT, p::FT, ρ::FT, ρe::FT)
vars_state(::DryAtmosModel, ::NoReferenceState, ::Auxiliary, FT) = @vars()

function init_state_auxiliary!(
    m::DryAtmosModel,
    ref_state::DryReferenceState,
    state_auxiliary,
    geom,
)
    FT = eltype(state_auxiliary)
    z = altitude(m, m.orientation, geom)
    T, p = ref_state.temperature_profile(param_set, z)

    _R_d::FT = R_d(param_set)
    ρ = p / (_R_d * T)
    Φ = state_auxiliary.Φ
    ρu = SVector{3, FT}(0, 0, 0)

    state_auxiliary.ref_state.T = T
    state_auxiliary.ref_state.p = p
    state_auxiliary.ref_state.ρ = ρ
    state_auxiliary.ref_state.ρe = totalenergy(ρ, ρu, p, Φ)
end

@inline function flux_first_order!(
    m::DryAtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    ρinv = 1 / ρ
    ρu = state.ρu
    ρe = state.ρe
    u = ρinv * ρu
    Φ = aux.Φ

    p = pressure(ρ, ρu, ρe, Φ)

    flux.ρ = ρ * u
    flux.ρu = p * I + ρ * u .* u'
    flux.ρe = u * (state.ρe + p)
end

function wavespeed(
    ::DryAtmosModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    p = pressure(ρ, ρu, ρe, Φ)

    u = ρu / ρ
    uN = abs(dot(nM, u))
    return uN + soundspeed(ρ, p)
end

"""
    pressure(ρ, ρu, ρe, Φ)

Compute the pressure given density `ρ`, momentum `ρu`, total energy `ρe`, and
gravitational potential `Φ`.
"""
function pressure(ρ, ρu, ρe, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ - ρ * Φ)
    else
        (γ - 1) * (ρe - dot(ρu, ρu) / 2ρ)
    end
end

"""
    totalenergy(ρ, ρu, p, Φ)

Compute the total energy given density `ρ`, momentum `ρu`, pressure `p`, and
gravitational potential `Φ`.
"""
function totalenergy(ρ, ρu, p, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        return p / (γ - 1) + dot(ρu, ρu) / 2ρ + ρ * Φ
    else
        return p / (γ - 1) + dot(ρu, ρu) / 2ρ
    end
end

"""
    soundspeed(ρ, p)

Compute the speed of sound from the density `ρ` and pressure `p`.
"""
function soundspeed(ρ, p)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    sqrt(γ * p / ρ)
end

"""
    vars_state(::DryAtmosModel, ::Prognostic, FT)

The prognostic state variables for the `DryAtmosModel` are density `ρ`, momentum `ρu`,
and total energy `ρe`
"""
function vars_state(::DryAtmosModel, ::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
    end
end

"""
    vars_state(::DryAtmosModel, ::Auxiliary, FT)

The auxiliary variables for the `DryAtmosModel` is gravitational potential
`Φ`
"""
function vars_state(m::DryAtmosModel, st::Auxiliary, FT)
    @vars begin
        Φ::FT
        ∇Φ::SVector{3, FT} # TODO: only needed for the linear model
        ref_state::vars_state(m, m.ref_state, st, FT)
        problem::vars_state(m, m.problem, st, FT)
    end
end
vars_state(::DryAtmosModel, ::AbstractDryAtmosProblem, ::Auxiliary, FT) =
    @vars()

"""
    vars_state(::DryAtmosModel, ::Entropy, FT)

The entropy variables for the `DryAtmosModel` correspond to the state
variables density `ρ`, momentum `ρu`, and total energy `ρe` as well as the
auxiliary variable gravitational potential `Φ`
"""
function vars_state(::DryAtmosModel, ::Entropy, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        Φ::FT
    end
end

"""
    state_to_entropy_variables!(
        ::DryAtmosModel,
        entropy::Vars,
        state::Vars,
        aux::Vars,
    )

See [`BalanceLaws.state_to_entropy_variables!`](@ref)
"""
function state_to_entropy_variables!(
    ::DryAtmosModel,
    entropy::Vars,
    state::Vars,
    aux::Vars,
)
    ρ, ρu, ρe, Φ = state.ρ, state.ρu, state.ρe, aux.Φ

    FT = eltype(state)
    γ = FT(gamma(param_set))

    p = pressure(ρ, ρu, ρe, Φ)
    s = log(p / ρ^γ)
    b = ρ / 2p
    u = ρu / ρ

    if total_energy
        entropy.ρ = (γ - s) / (γ - 1) - (dot(u, u) - 2Φ) * b
    else
        entropy.ρ = (γ - s) / (γ - 1) - (dot(u, u)) * b
    end
    entropy.ρu = 2b * u
    entropy.ρe = -2b
    entropy.Φ = 2ρ * b
end

"""
    entropy_variables_to_state!(
        ::DryAtmosModel,
        state::Vars,
        aux::Vars,
        entropy::Vars,
    )

See [`BalanceLaws.entropy_variables_to_state!`](@ref)
"""
function entropy_variables_to_state!(
    ::DryAtmosModel,
    state::Vars,
    aux::Vars,
    entropy::Vars,
)
    FT = eltype(state)
    β = entropy
    γ = FT(gamma(param_set))

    b = -β.ρe / 2
    ρ = β.Φ / (2b)
    ρu = ρ * β.ρu / (2b)

    p = ρ / (2b)
    s = log(p / ρ^γ)
    Φ = dot(ρu, ρu) / (2 * ρ^2) - ((γ - s) / (γ - 1) - β.ρ) / (2b)

    ρe = p / (γ - 1) + dot(ρu, ρu) / (2ρ) + ρ * Φ

    state.ρ = ρ
    state.ρu = ρu
    state.ρe = ρe
    aux.Φ = Φ
end

function state_to_entropy(::DryAtmosModel, state::Vars, aux::Vars)
    FT = eltype(state)
    ρ, ρu, ρe, Φ = state.ρ, state.ρu, state.ρe, aux.Φ
    p = pressure(ρ, ρu, ρe, Φ)
    γ = FT(gamma(param_set))
    s = log(p / ρ^γ)
    η = -ρ * s / (γ - 1)
    return η
end

function numerical_volume_conservative_flux_first_order!(
    ::EntropyConservative,
    ::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    FT = eltype(F)
    ρ_1, ρu_1, ρe_1 = state_1.ρ, state_1.ρu, state_1.ρe
    ρ_2, ρu_2, ρe_2 = state_2.ρ, state_2.ρu, state_2.ρe
    Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ
    u_1 = ρu_1 / ρ_1
    u_2 = ρu_2 / ρ_2
    p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)
    p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)
    b_1 = ρ_1 / 2p_1
    b_2 = ρ_2 / 2p_2

    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    b_avg = ave(b_1, b_2)
    Φ_avg = ave(Φ_1, Φ_2)

    usq_avg = ave(dot(u_1, u_1), dot(u_2, u_2))

    ρ_log = logave(ρ_1, ρ_2)
    b_log = logave(b_1, b_2)
    α = b_avg * ρ_log / 2b_1

    γ = FT(gamma(param_set))

    Fρ = u_avg * ρ_log
    Fρu = u_avg * Fρ' + ρ_avg / 2b_avg * I
    if total_energy
        Fρe =
            (1 / (2 * (γ - 1) * b_log) - usq_avg / 2 + Φ_avg) * Fρ + Fρu * u_avg
    else
        Fρe = (1 / (2 * (γ - 1) * b_log) - usq_avg / 2) * Fρ + Fρu * u_avg
    end

    F.ρ += Fρ
    F.ρu += Fρu
    F.ρe += Fρe
end

function numerical_volume_fluctuation_flux_first_order!(
    ::NumericalFluxFirstOrder,
    ::DryAtmosModel,
    D::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    if fluctuation_gravity
        FT = eltype(D)
        ρ_1, ρu_1, ρe_1 = state_1.ρ, state_1.ρu, state_1.ρe
        ρ_2, ρu_2, ρe_2 = state_2.ρ, state_2.ρu, state_2.ρe
        Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ
        p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)
        p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)
        b_1 = ρ_1 / 2p_1
        b_2 = ρ_2 / 2p_2

        ρ_log = logave(ρ_1, ρ_2)
        b_avg = ave(b_1, b_2)
        α = b_avg * ρ_log / 2b_1

        D.ρu -= α * (Φ_1 - Φ_2) * I
    end
end

struct CentralVolumeFlux <: NumericalFluxFirstOrder end
function numerical_volume_conservative_flux_first_order!(
    ::CentralVolumeFlux,
    m::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    FT = eltype(F)
    F_1 = similar(F)
    flux_first_order!(m, F_1, state_1, aux_1, FT(0), EveryDirection())

    F_2 = similar(F)
    flux_first_order!(m, F_2, state_2, aux_2, FT(0), EveryDirection())

    parent(F) .= (parent(F_1) .+ parent(F_2)) ./ 2
end

struct KGVolumeFlux <: NumericalFluxFirstOrder end
function numerical_volume_conservative_flux_first_order!(
    ::KGVolumeFlux,
    m::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    Φ_1 = aux_1.Φ
    ρ_1 = state_1.ρ
    ρu_1 = state_1.ρu
    ρe_1 = state_1.ρe
    u_1 = ρu_1 / ρ_1
    e_1 = ρe_1 / ρ_1
    p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)

    Φ_2 = aux_2.Φ
    ρ_2 = state_2.ρ
    ρu_2 = state_2.ρu
    ρe_2 = state_2.ρe
    u_2 = ρu_2 / ρ_2
    e_2 = ρe_2 / ρ_2
    p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)

    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    e_avg = ave(e_1, e_2)
    p_avg = ave(p_1, p_2)

    F.ρ = ρ_avg * u_avg
    F.ρu = p_avg * I + ρ_avg * u_avg .* u_avg'
    F.ρe = ρ_avg * u_avg * e_avg + p_avg * u_avg
end

struct Coriolis end
function source!(
    m::DryAtmosModel,
    ::Coriolis,
    source,
    state_prognostic,
    state_auxiliary,
)
    FT = eltype(state_prognostic)
    _Omega::FT = Omega(param_set)
    # note: this assumes a SphericalOrientation
    source.ρu -= SVector(0, 0, 2 * _Omega) × state_prognostic.ρu
end

function source!(m::DryAtmosModel, source, state_prognostic, state_auxiliary)
    ntuple(Val(length(m.sources))) do s
        Base.@_inline_meta
        source!(m, m.sources[s], source, state_prognostic, state_auxiliary)
    end
end

struct EntropyConservativeWithPenalty <: NumericalFluxFirstOrder end
function numerical_flux_first_order!(
    numerical_flux::EntropyConservativeWithPenalty,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    FT = eltype(fluxᵀn)

    numerical_flux_first_order!(
        EntropyConservative(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    fluxᵀn = parent(fluxᵀn)

    wavespeed⁻ = wavespeed(
        balance_law,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        direction,
    )
    wavespeed⁺ = wavespeed(
        balance_law,
        normal_vector,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    max_wavespeed = max.(wavespeed⁻, wavespeed⁺)
    penalty =
        max_wavespeed .* (parent(state_prognostic⁻) - parent(state_prognostic⁺))

    fluxᵀn .+= penalty / 2
end

Base.@kwdef struct MatrixFlux{FT} <: NumericalFluxFirstOrder
    Mcut::FT = 0
    low_mach::Bool = false
    kinetic_energy_preserving::Bool = false
end

function numerical_flux_first_order!(
    numerical_flux::MatrixFlux,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    FT = eltype(fluxᵀn)
    numerical_flux_first_order!(
        EntropyConservative(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    fluxᵀn = parent(fluxᵀn)

    γ = FT(gamma(param_set))

    low_mach = numerical_flux.low_mach
    Mcut = numerical_flux.Mcut
    kinetic_energy_preserving = numerical_flux.kinetic_energy_preserving

    ω = FT(π) / 3
    δ = FT(π) / 5
    random_unit_vector = SVector(sin(ω) * cos(δ), cos(ω) * cos(δ), sin(δ))
    # tangent space basis
    τ1 = random_unit_vector × normal_vector
    τ2 = τ1 × normal_vector

    ρ⁻ = state_prognostic⁻.ρ
    ρu⁻ = state_prognostic⁻.ρu
    ρe⁻ = state_prognostic⁻.ρe

    Φ⁻ = state_auxiliary⁻.Φ
    u⁻ = ρu⁻ / ρ⁻
    p⁻ = pressure(ρ⁻, ρu⁻, ρe⁻, Φ⁻)
    β⁻ = ρ⁻ / 2p⁻

    Φ⁺ = state_auxiliary⁺.Φ
    ρ⁺ = state_prognostic⁺.ρ
    ρu⁺ = state_prognostic⁺.ρu
    ρe⁺ = state_prognostic⁺.ρe

    u⁺ = ρu⁺ / ρ⁺
    p⁺ = pressure(ρ⁺, ρu⁺, ρe⁺, Φ⁺)
    β⁺ = ρ⁺ / 2p⁺

    ρ_log = logave(ρ⁻, ρ⁺)
    β_log = logave(β⁻, β⁺)
    if total_energy
        Φ_avg = ave(Φ⁻, Φ⁺)
    else
        Φ_avg = 0
    end
    u_avg = ave.(u⁻, u⁺)
    p_avg = ave(ρ⁻, ρ⁺) / 2ave(β⁻, β⁺)
    u²_bar = 2 * sum(u_avg .^ 2) - sum(ave(u⁻ .^ 2, u⁺ .^ 2))

    h_bar = γ / (2 * β_log * (γ - 1)) + u²_bar / 2 + Φ_avg
    c_bar = sqrt(γ * p_avg / ρ_log)

    umc = u_avg - c_bar * normal_vector
    upc = u_avg + c_bar * normal_vector
    u_avgᵀn = u_avg' * normal_vector
    R = hcat(
        SVector(1, umc[1], umc[2], umc[3], h_bar - c_bar * u_avgᵀn),
        SVector(1, u_avg[1], u_avg[2], u_avg[3], u²_bar / 2 + Φ_avg),
        SVector(0, τ1[1], τ1[2], τ1[3], τ1' * u_avg),
        SVector(0, τ2[1], τ2[2], τ2[3], τ2' * u_avg),
        SVector(1, upc[1], upc[2], upc[3], h_bar + c_bar * u_avgᵀn),
    )

    if low_mach
        M = abs(u_avg' * normal_vector) / c_bar
        c_bar *= max(min(M, FT(1)), Mcut)
    end

    if kinetic_energy_preserving
        λl = abs(u_avgᵀn) + c_bar
        λr = λl
    else
        λl = abs(u_avgᵀn - c_bar)
        λr = abs(u_avgᵀn + c_bar)
    end

    Λ = SDiagonal(λl, abs(u_avgᵀn), abs(u_avgᵀn), abs(u_avgᵀn), λr)
    #Ξ = sqrt(abs((p⁺ - p⁻) / (p⁺ + p⁻)))
    #Λ = Ξ * abs(u_avgᵀn + c_bar) * I + (1 - Ξ) * ΛM

    T = SDiagonal(ρ_log / 2γ, ρ_log * (γ - 1) / γ, p_avg, p_avg, ρ_log / 2γ)

    entropy⁻ = similar(parent(state_prognostic⁻), Size(6))
    state_to_entropy_variables!(
        balance_law,
        Vars{vars_state(balance_law, Entropy(), FT)}(entropy⁻),
        state_prognostic⁻,
        state_auxiliary⁻,
    )

    entropy⁺ = similar(parent(state_prognostic⁺), Size(6))
    state_to_entropy_variables!(
        balance_law,
        Vars{vars_state(balance_law, Entropy(), FT)}(entropy⁺),
        state_prognostic⁺,
        state_auxiliary⁺,
    )

    Δentropy = parent(entropy⁺) - parent(entropy⁻)

    fluxᵀn .-= R * Λ * T * R' * Δentropy[SOneTo(5)] / 2
end

function vertical_unit_vector(m::DryAtmosModel, aux::Vars)
    FT = eltype(aux)
    aux.∇Φ / FT(grav(param_set))
end

norm_u(state::Vars, k̂::AbstractVector, ::VerticalDirection) =
    abs(dot(state.ρu, k̂)) / state.ρ
norm_u(state::Vars, k̂::AbstractVector, ::HorizontalDirection) =
    norm((state.ρu .- dot(state.ρu, k̂) * k̂) / state.ρ)
norm_u(state::Vars, k̂::AbstractVector, ::Direction) = norm(state.ρu / state.ρ)

function advective_courant(
    m::DryAtmosModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    k̂ = vertical_unit_vector(m, aux)
    normu = norm_u(state, k̂, direction)
    return Δt * normu / Δx
end

function nondiffusive_courant(
    m::DryAtmosModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    p = pressure(ρ, ρu, ρe, Φ)

    k̂ = vertical_unit_vector(m, aux)
    normu = norm_u(state, k̂, direction)
    ss = soundspeed(ρ, p)
    return Δt * (normu + ss) / Δx
end

function drag_source!(m::DryAtmosModel, args...)
    drag_source!(m, m.drag_source, args...)
end
struct NoDrag end
drag_source!(m::DryAtmosModel, ::NoDrag, args...) = nothing

struct Gravity end
function source!(m::DryAtmosModel, ::Gravity, source, state, aux)
    ∇Φ = aux.∇Φ
    if !fluctuation_gravity
        source.ρu -= state.ρ * ∇Φ
    end
    if !total_energy
        source.ρe -= state.ρu' * ∇Φ
    end
end

# Numerical Flux
#=
numerical_flux_first_order!(::Nothing, _...) = nothing
numerical_boundary_flux_first_order!(::Nothing, _...) = nothing
numerical_flux_second_order!(::Nothing, _...) = nothing
numerical_boundary_flux_second_order!(::Nothing, _...) = nothing
=#
numerical_flux_first_order!(::Nothing, ::DryAtmosModel, _...) = nothing
numerical_flux_second_order!(::Nothing, ::DryAtmosModel, _...) = nothing
numerical_boundary_flux_second_order!(::Nothing, a, ::DryAtmosModel, _...) =
    nothing

include("linear.jl")

# Throwing this in for convenience
function cubedshellwarp(a, b, c, R = max(abs(a), abs(b), abs(c)))

    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        x1 = sR / sqrt(X^2 + Y^2 + 1)
        x2, x3 = X * x1, Y * x1
        x1, x2, x3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : Face I from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : Face II from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : Face III from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : Face IV from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : Face V from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : Face VI from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubedshellwarp: $a, $b, $c")
    end

    return x1, x2, x3
end


text:
export ResNN,getResNN

"""
Residual Neural Network block

Y_k+1 = Y_k + h*layer{k}(theta{k},Y_k)
"""
mutable struct ResNN{T, TL <: AbstractMeganetElement{T}} <: AbstractMeganetElement{T} #TODO limit TL more
    layer    :: TL
    nt       :: Int
    h        :: T
    outTimes :: Array{Int,1}
    Q
end

function getResNN(TYPE::Type,layer,nt,h=one(TYPE),outTimes=eye(Int,nt)[:,nt],Q=I)
    h = convert(TYPE,h);
	if nFeatIn(layer)!=nFeatOut(layer)
            error("ResNN layer must be square!")
    end
    return ResNN(layer,nt,h,outTimes,Q)
end

function nTheta(this::ResNN)
    return this.nt*nTheta(this.layer);
end

function nFeatIn(this::ResNN)
    return nFeatIn(this.layer);
end

function nFeatOut(this::ResNN)
    return nFeatOut(this.layer);
end

function nDataOut(this::ResNN)
    if length(this.Q)==1
        n = sum(this.outTimes)*nFeatOut(this.layer)
    else
        n = sum(this.outTimes)*size(this.Q,1)
    end
    return n
end

function initTheta(this::ResNN{T}) where {T<:Number}
    return repmat(vec(initTheta(this.layer)),this.nt,1)
end

# ------- apply forward problems -----------
function  apply(this::ResNN{T},theta_in::Array{T},Y0::Array{T},tmp,doDerivative=true) where {T<:Number}
    if isempty(tmp)
        tmp = Array{Any}(this.nt+1,2)
    end

    if doDerivative
        if isassigned(tmp,1,1)
            tmp11 = tmp[1,1]
            tmp11 .= Y0
        else
            tmp[1,1] = copy(Y0)
        end
    end

    nex = div(length(Y0),nFeatIn(this))
    Y   = reshape(Y0,:,nex)

    theta = reshape(theta_in,:,this.nt)

    Ydata::Array{T,2} = zeros(T,0,nex)
    for i=1:this.nt
        if !isassigned(tmp,i,2)
            tmp[i,2] = Array{Any}(0)
        end
        Z,dummy,tmp[i,2] = apply(this.layer,theta[:,i],Y,tmp[i,2],doDerivative)
        Y +=  this.h * Z
        if this.outTimes[i]==1
            Ydata = [Ydata;this.Q*Y]
        end

        if doDerivative
            if isassigned(tmp,i+1,1)
                tmp1 = tmp[i+1,1]
                tmp1 .= Y
            else
                tmp[i+1,1] = copy(Y)
            end
        end

    end
    return Ydata,Y,tmp
end

# -------- Jacobian matvecs ---------------
function JYmv(this::ResNN{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T<:Number}
    # if isempty(dY)
    #     dY = 0.0;
    # elseif length(dY)>1
        nex = div(length(dY),nFeatIn(this))
        dY   = reshape(dY,:,nex)
    # end
    dYdata = zeros(T,0,nex);
    theta  = reshape(theta,:,this.nt);
    for i=1:this.nt
        dY += this.h* JYmv(this.layer,dY,theta[:,i],tmp[i,1],tmp[i,2])[2];
        if this.outTimes[i]==1
            dYdata = [dYdata; this.Q*dY];
        end
    end
    return dYdata, dY
end

function  Jmv(this::ResNN{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T<:Number}
    nex = div(length(Y),nFeatIn(this))
    if length(dY)==0
         dY = zeros(T,size(Y))
     elseif length(dY)>1
        dY   = reshape(dY,:,nex)
    end

    dYdata = zeros(T,0,nex)
    theta  = reshape(theta,:,this.nt)
    dtheta = reshape(dtheta,:,this.nt)
    for i=1:this.nt
        dY +=  this.h* Jmv(this.layer,dtheta[:,i],dY,theta[:,i],tmp[i,1],tmp[i,2])[2]
        if this.outTimes[i]==1
            dYdata = [dYdata;this.Q*dY]
        end
    end
    return dYdata,dY
end

# -------- Jacobian transpose matvecs ----------------

function JYTmv(this::ResNN{T},Wdata_in::Array{T},Win::Array{T},theta_in::Array{T},Y::Array{T},tmp) where {T<:Number}

    nex = div(length(Y),nFeatIn(this))
    if length(Wdata_in)>0
        Wdata = reshape(Wdata_in,:,sum(this.outTimes),nex)
    end
    if length(Win)==0
        W = zeros(T,0,0)
    else
        W = reshape(Win,:,nex)
    end
    theta  = reshape(theta_in,:,this.nt)

    cnt = sum(this.outTimes)
    for i=this.nt:-1:1
        if  this.outTimes[i]==1
            W::Array{T,2} = this.Q'*Wdata[:,cnt,:]
            cnt = cnt-1
        end
        dW::Array{T,2} = JYTmv(this.layer,W,zeros(T,0),theta[:,i],tmp[i,1],tmp[i,2])
        W  += this.h*dW
    end
    return W
end

function JTmv(this::ResNN{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T<:Number}

    nex = div(length(Y),nFeatIn(this))
    if !isempty(Wdata) && any(this.outTimes.!=0)
        Wdata = reshape(Wdata,:,sum(this.outTimes),nex)
    end
    if length(W)==0
        if any(this.outTimes.!=0)
            W = zero(T) #assume the Wdata is non-zero
        else
            W = Wdata
        end
    end
    if length(W)>1
        W     = reshape(W,:,nex)
    end

    theta  = reshape(theta,:,this.nt)
    cnt    = sum(this.outTimes)
    dtheta = zeros(T,size(theta))

    for i=this.nt:-1:1
        if  this.outTimes[i]==1
            W +=  this.Q'* Wdata[:,cnt,:]
            cnt = cnt-1
        end

        dmbi,dW = JTmv(this.layer,W,zeros(T,0),theta[:,i],tmp[i,1],tmp[i,2])

        dtheta[:,i]  = this.h*dmbi
        W += this.h*dW
    end

    return vec(dtheta), W
end


text:
a = 42

@doc raw"""
    func2big(x)

Computes the square of `x` and add `a` to it.

# Examples
```julia-repl
julia> include("big-module.jl")
julia> using .MyBigModule
julia> func2big(3)
51
```
"""
function func2big(x)
    return func1big(x) + a
end

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
abstract type OptimizationContainerKey end

const _DELIMITER = "__"

function make_key(::Type{T}, args...) where {T <: OptimizationContainerKey}
    return T(args...)
end

function encode_key(key::OptimizationContainerKey)
    return encode_symbol(get_component_type(key), get_entry_type(key), key.meta)
end

encode_key_as_string(key::OptimizationContainerKey) = string(encode_key(key))
encode_keys_as_strings(container_keys) = [encode_key_as_string(k) for k in container_keys]

function encode_symbol(
    ::Type{T},
    ::Type{U},
    meta::String=CONTAINER_KEY_EMPTY_META,
) where {T <: Union{PSY.Component, PSY.System}, U}
    meta_ = isempty(meta) ? meta : _DELIMITER * meta
    T_ = replace(replace(IS.strip_module_name(T), "{" => _DELIMITER), "}" => "")
    return Symbol("$(IS.strip_module_name(string(U)))$(_DELIMITER)$(T_)" * meta_)
end

function check_meta_chars(meta)
    # Underscores in this field will prevent us from being able to decode keys.
    if occursin(_DELIMITER, meta)
        throw(IS.InvalidValue("'$_DELIMITER' is not allowed in meta"))
    end
end

function should_write_resulting_value(key_val::OptimizationContainerKey)
    value_type = get_entry_type(key_val)
    return should_write_resulting_value(value_type)
end

function convert_result_to_natural_units(key::OptimizationContainerKey)
    return convert_result_to_natural_units(get_entry_type(key))
end


text:
##
"""
CELL5M is a 5-minute arc resolution view of the -180,-90
lat-long system. This corresponds to a cell-size of
0.08333333333. The following code converts between cell5m
values and latlong.
"""
function cell5m_to_latlong(cell::Int)
  nrows = 4320
  ncols = 2160
  col, row = divrem(cell, nrows)
  x = row*0.0833333333333333333-(180-0.0833333333333333333/2)
  y = (90-0.0833333333333333333/2) - col*0.0833333333333333333
  return x,y
end
cell5m_to_latlong(1300706)
## Test data (from spam2010V1r1_global_H_TA)

cell5m_test = [1300706
1300707
1300708
1305026
1305027
1305028
1348258
1348259
1348260
1352579
1352580
1473490]
cell5m_ref_x = [-147.7916666668
-147.7083333335
-147.6250000001
-147.7916666668
-147.7083333335
-147.6250000001
-145.1250000001
-145.0416666668
-144.9583333335
-145.0416666668
-144.9583333335
-149.1250000001]
cell5m_ref_y = [ 64.8749999994
 64.8749999994
 64.8749999994
 64.791666666
 64.791666666
 64.791666666
 63.9583333327
 63.9583333327
 63.9583333327
 63.8749999994
 63.8749999994
 61.5416666661]
for i in 1:length(cell5m_test)
  x,y = cell5m_to_latlong(cell5m_test[i])
  println(cell5m_test[i], " ",
          x, " ", cell5m_ref_x[i], " ",
          y, " ", cell5m_ref_y[i])
  @assert (x ≈ cell5m_ref_x[i])
  @assert (y ≈ cell5m_ref_y[i])
end
##

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
struct PValue{TI<:Integer,TR<:Real,TO<:Ordering,TD<:AbstractDesign}
    ordering::TO
    design::TD
    XX::Array{TI,2}
    inds::Array{Bool,2}
    p0::TR
end

Base.iterate(pval::PValue, state = 0) = state > 0 ? nothing : (pval, state + 1)
Base.length(pval::PValue) = 1

Base.show(io::IO, pval::PValue) = print(io, string(pval))
Base.show(io::IO, ::MIME"application/prs.juno.inline", pval::PValue) = print(io, string(pval))
string(pval::PValue) = "p-value"

function PValue(ordering::Ordering, design::AbstractDesign, p0::Real)

    XX   = sample_space(design)
    nn   = size(XX, 1)
    inds = zeros(nn, nn)
    for i = 1:nn
        x1, x2    = XX[i,1], XX[i,2]
        inds[i,:] = more_extreme.(XX[:,1], XX[:,2], x1, x2, ordering, design)
    end
    return PValue{eltype(XX),typeof(p0),typeof(ordering),typeof(design)}(ordering, design, XX, inds, p0)
end

function PValue(estimator::Estimator, design::AbstractDesign, p0::Real; orientation::Symbol = :superiority)
    PValue(EstimatorOrdering(estimator), design, p0)
end

function (pval::PValue{TI,TR,TO,TD})(x1::TI, x2::TI) where {TI<:Integer,TR<:Real,TO<:Ordering,TD<:AbstractDesign}

    design = pval.design
    p0     = pval.p0
    XX     = pval.XX
    ind    = mapslices(all, XX .== [x1 x2], dims = 2) |> vec |> findall
    @assert length(ind) == 1
    ind    = ind[1]
    inds   = pval.inds[ind,:]
    return min(1, max(0, sum(
        pmf.(XX[inds,2], n2.(design, XX[inds,1]), p0) .*
        pmf.(XX[inds,1], n1(design), p0)
    ) ) )
end

evaluate(pval::PValue{TI,TR,TO,TD}, x1::TI, x2::TI) where {TI<:Integer,TR<:Real,TO<:Ordering,TD<:AbstractDesign} = pval(x1, x2)


text:
# Check for Model equality up to reorderings of a few fields
function ≊(m1::Model,m2::Model)
    function eq_tuples(nt1::NamedTuple,nt2::NamedTuple)
        return length(nt1)==length(nt2) && all(nt1[k]==nt2[k] for k in keys(nt1))
    end
    return Set(arguments(m1))==Set(arguments(m2)) && m1.retn==m2.retn && eq_tuples(m1.dists,m2.dists) && eq_tuples(m1.vals,m2.vals)
end

m = @model (n,α,β) begin
    p ~ Beta(α, β)
    x ~ Binomial(n, p)
    z ~ Binomial(n, α/(α+β))
end

@testset "prior" begin
    m1 = Soss.prior(m, :x)
    @test Soss.prior(m, :x) ≊ @model (α, β) begin
        p ~ Beta(α, β)
    end
end

@testset "likelihood" begin
    m1 = Soss.likelihood(m, :x)
    @test Soss.likelihood(m, :x) ≊ @model (p, n) begin
        x ~ Binomial(n, p)
    end
end

m1 = prune(m, :z)
@testset "prune" begin
    @test prune(m, :x, :z) ≊ @model (α, β) begin
        p ~ Beta(α, β)
    end
    @test prune(m1, :n) ≊ @model (α, β) begin
        p ~ Beta(α, β)
    end
    @test prune(m, :p) ≊ @model (α, n, β) begin
        z ~ Binomial(n, α / (α + β))
    end

    # When I define these variables, the tests pass.
    # Doing "@test prune(m1, :p) ≊ @model begin end" strangely causes an error in @model about reducing over an empty collection.
    emptymodel = @model begin end
    @test prune(m1, :p) ≊ emptymodel
end

@testset "predictive" begin
    @test predictive(m, :p) ≊ @model (n, p) begin
        x ~ Binomial(n, p)
    end
end

@testset "Do" begin
    @test Do(m, :p, :z) ≊ @model (n, p) begin
        x ~ Binomial(n, p)
    end
    empty = @model begin end
    @test Do(m, variables(m)...) ≊ empty
end


text:
import CLBLAS

const clblas = CLBLAS
clblas.setup()

using OpenCL: cl

device, ctx, queue = clblas.get_next_compute_context()

N = Unsigned(5)
data = Vector{Float32}(5)

data[1] = 5.55555
data[2] = 3.333
data[3] = 2.22
data[4] = 4.4444
data[5] = 1.1

bufX = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=data)
scratchBuff = cl.Buffer(Float32, ctx, :rw, 2 * length(data))
iMax = cl.Buffer(cl.CL_uint, ctx, :w, 1)

offiMax = offx = Unsigned(0)
incx = cl.cl_int(1)
ncq = cl.cl_uint(1)
clq = queue.id
newl = cl.cl_uint(0)

event = cl.UserEvent(ctx, retain=true)
ptrEvent = [event.id]

clblas.clblasiSamax( N, iMax.id, offiMax, bufX.id, offx, incx, scratchBuff.id, ncq, [clq], newl, C_NULL, ptrEvent);

cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

result = Array(cl.CL_uint, 1)
cl.enqueue_read_buffer(queue, iMax, result, Unsigned(0), nothing, true)

result = int(floor(result[1]))

expected = float32(5.55555)
println(data[result])
println("------------")
println(expected)

try
    @assert(isapprox(norm(expected - data[result]), zero(Float32)))
    info("success!")
finally
    clblas.teardown()
end


text:
compress_then_write(b::CategoricalVector{T, IntType}, io) where {T, IntType <: Integer} = begin
    #println("abc")
    compress_refs = compress_then_write(b.refs, io)
    compress_poolindex = compress_then_write(b.pool.index, io)

    (type = CategoricalVector, refs = compress_refs, poolindex = compress_poolindex, ordered = b.pool.ordered)
end

column_loader(b::Type{CategoricalVector}, io, metadata) = begin
    refs_meta = metadata.refs
    pi_meta = metadata.poolindex
    ref = column_loader(refs_meta.type, io, refs_meta)
    poolindex = column_loader(pi_meta.type, io, pi_meta)
    CategoricalArray{pi_meta.type, 1}(ref, CategoricalPool{eltype(poolindex), eltype(ref)}(Array(poolindex), metadata.ordered))
end

if false
# # compress_then_write-cate.jl
# using Revise
# using JDF
# using DataFrames, Debugger
# b = categorical(["a", "b", "a", "c"])
# io = open("c:/data/test.io", "w")
# metadata = compress_then_write(a, io)
# close(io)
#
# io = open("c:/data/test.io", "r")
# aa = column_loader(CategoricalVector, io, metadata)
# close(io)
#
# df = DataFrame(a = a)
#
# savejdf("c:/data/pls_del.jdf", df)
#
# loadjdf("c:/data/pls_del.jdf")
end


text:
module Nanosoldier

using Dates, Distributed, Printf, InteractiveUtils
import GitHub, BenchmarkTools, JSON, HTTP

const SHA_SEPARATOR = '@'
const BRANCH_SEPARATOR = ':'
const TAG_SEPARATOR = '#'

#####################
# utility functions #
#####################

snip(str, len) = str[1:min(len, end)]
snipsha(sha) = snip(sha, 7)

function gitclone!(repo, path, auth=nothing)
    if isa(auth, GitHub.OAuth2)
        run(`git clone https://$(auth.token):x-oauth-basic@github.com/$(repo).git $(path)`)
    elseif isa(auth, GitHub.UsernamePassAuth)
        run(`git clone https://$(auth.username):$(auth.password)@github.com/$(repo).git $(path)`)
    else
        run(`git clone git@github.com:$(repo).git $(path)`)
    end
end

gitreset!() = (run(`git fetch --all`); run(`git reset --hard origin/master`))
gitreset!(path) = cd(gitreset!, path)

##################
# error handling #
##################

mutable struct NanosoldierError{E<:Exception} <: Exception
    url::String
    msg::String
    err::E
end

NanosoldierError(msg, err::E) where {E<:Exception} = NanosoldierError{E}("", msg, err)

function Base.show(io::IO, err::NanosoldierError)
    print(io, "NanosoldierError: ", err.msg, ": ")
    showerror(io, err.err)
end

############
# includes #
############

include("config.jl")
include("build.jl")
include("submission.jl")
include("jobs/jobs.jl")
include("server.jl")

end # module


text:
export TR_Sec_ls
"""
A one dimensionnal trust region algorithm using the secant as the approximation
for the second derivative in the second order Taylor expansion.
For more documentation referer to TR_generic_ls.
"""
function TR_Sec_ls(h :: LineModel,  stop_ls :: LS_Stopping;  kwargs...)
    ls_meta = LS_Function_Meta(dir = :Sec)
    (state, stop_ls.meta.optimal) = TR_generic_ls(h, stop_ls; f_meta = ls_meta, kwargs...)
    return (state, stop_ls.meta.optimal)
end


text:
using ContinuumArrays, Test
import ContinuumArrays: AffineQuasiVector

@testset "AffineQuasiVector" begin
    x = Inclusion(0..1)
    y = 2x .- 1
    a = affine(-2..3, -1..1)

    @testset "algebra" begin
        @test 2x isa AffineQuasiVector
        @test x/2 isa AffineQuasiVector
        @test (2x)[0.1] == 0.2
        @test_throws BoundsError (2x)[2]

        @test y isa AffineQuasiVector
        @test y[0.1] == 2*(0.1)-1
        @test y/2 isa AffineQuasiVector
        @test 2\y isa AffineQuasiVector
        @test (y/2)[0.1] == (2\y)[0.1] == -0.4
        @test y .+ 1 isa AffineQuasiVector
        @test (y .+ 1)[0.1] == (1 .+ y)[0.1]
        @test (y .- 1)[0.1] == y[0.1]-1
        @test (1 .- y)[0.1] == 1-y[0.1]
        @test y[x] == y[:] == y
        @test_throws BoundsError y[Inclusion(0..2)]

        @test y == y

        @test x == affine(x,x)
        @test affine(x,x) == x
    end

    @testset "find" begin
        @test findfirst(isequal(0.1),x) == findlast(isequal(0.1),x) == 0.1
        @test findall(isequal(0.1), x) == [0.1]
        @test findfirst(isequal(2),x) == findlast(isequal(2),x) == nothing
        @test findall(isequal(2), x) == Float64[]

        @test findfirst(isequal(0.2),y) == findlast(isequal(0.2),y) == 0.6
        @test findfirst(isequal(2.3),y) == findlast(isequal(2.3),y) == nothing
        @test findall(isequal(0.2),y) == [0.6]
        @test findall(isequal(2),y) == Float64[]
    end

    @testset "minmax" begin
        @test AffineQuasiVector(x)[0.1] == 0.1
        @test minimum(y) == -1
        @test maximum(y) == 1
        @test union(y) == Inclusion(-1..1)
        @test ContinuumArrays.inbounds_getindex(y,0.1) == y[0.1]
        @test ContinuumArrays.inbounds_getindex(y,2.1) == 2*2.1 - 1

        z = 1 .- x
        @test minimum(z) == 0.0
        @test maximum(z) == 1.0
        @test union(z) == Inclusion(0..1)

        @test !isempty(z)
        @test z == z
    end

    @testset "AffineMap" begin
        @test a[0.1] == -0.16
        @test_throws BoundsError a[-3]
        @test a[-2] == first(a) == -1
        @test a[3] == last(a) == 1
        @test invmap(a)[-0.16] ≈ 0.1
        @test invmap(a)[1] == 3
        @test invmap(a)[-1] == -2
        @test union(a) == Inclusion(-1..1)

        @test affine(0..1, -1..1) == y
    end

    @testset "show" begin
        @test stringmime("text/plain", y) == "2.0 * Inclusion(0..1) .+ (-1.0)"
        @test stringmime("text/plain", a) == "0.4 * Inclusion(-2..3) .+ (-0.2)"
    end
end