text:
using Test
using UnixIO
using UnixIO: dbstring, C
using AsyncLog
using LoggingTestSets
using Crayons

@testset LoggingTestSet "Pseudoterminal Tests" begin

function ptdump(f, cmd, cin, cout)

    blue = crayon"fg:blue"
    red = crayon"bold fg:red"

    println("┌ ", blue, cmd, inv(blue), ":")
    pad =   "│ "
    needpad = Ref(false)
    neednl = Ref(false)

    function fin()
        x = ""
        for _ in 1:30
            x = readline(cout; keep=true, timeout=0.1)
            if x != ""
                break
            end
        end
        if x == ""
            return nothing
        end
        if neednl[]
            print(red, "!", inv(red), "\n")
        end
        neednl[] = false
        print(pad)
        print(blue, x, inv(blue))
        needpad[] = x[end] == '\n'
        x
    end

    function fout(x)
        if needpad[]
            print(pad)
        end
        for c in x
            write(cin, c)
            if c isa Char
                print(c)
            end
            sleep(0.05 * rand())
        end
        needpad[] = x[end] == '\n'
        neednl[] = x[end] != '\n'
    end

    f(fin, fout)

    println("\n└")
end

opts = (check_status=false,)

@info "Testing with Julia UnixIO pty clinet"

mktempdir() do d
    jlfile=joinpath(d, "tmp.jl")
    write(jlfile, """
        using UnixIO
        io = UnixIO.stdin
        @show ENV["JULIA_UNIX_IO_DEBUG_LEVEL"]
        @show UnixIO.DEBUG_LEVEL
        @show typeof(io)
        @assert io isa UnixIO.FD{UnixIO.In,UnixIO.CanonicalMode}
        @assert !(io isa UnixIO.FD{UnixIO.In,UnixIO.Pseudoterminal})
        @show UnixIO.ReadFragmentation(io)
        println("go!")
        while true
            l = readline(io; keep=true)
            b = bytesavailable(io.buffer)
            UnixIO.println(length(l), ":", repr(l), ":", b)
        end
        """)

    @info "readline() from pseudoterminal with fragmented writes."
    env = merge(ENV, Dict("JULIA_UNIX_IO_DEBUG_LEVEL" => "0",
                          "JULIA_DEBUG" => ""))
    UnixIO.ptopen(`$(Base.julia_cmd()) $jlfile`; env=env, opts...) do cin, cout
    ptdump(Base.julia_cmd(), cin, cout) do fin, fout
        while true
            l = fin()
            if l == nothing
                sleep(1)
            elseif contains(l, r"go!")
                break
            end
        end

        fout("Hello1\n")
        fout("Hello2\n")
        @test fin() == """7:"Hello1\\n":0\n"""
        @test fin() == """7:"Hello2\\n":0\n"""

        fout("Hel")

        # Even after draining and waiting, the client read() has not returned
        # because of canonical mode.
        @test fin() == nothing
        fout("lo3\n")
        @test fin() == """7:"Hello3\\n":0\n"""

        # `CEOF` and sends a line with no "\n" (e.g. like the "bash$ " prompt)
        fout("Hello4")
        fout([UInt8(C.CEOF)])
        @test fin() == """6:"Hello4":0\n"""
        fout("Hello5\x04Hello6\n")
        @test fin() == """6:"Hello5":0\n"""
        @test fin() == """7:"Hello6\\n":0\n"""

        fout("Hello7\r")
        @test fin() == nothing
        fout("\x04")
        @test fin() == """7:"Hello7\\r":0\n"""
    end
    end
end

@info "Testing with C blocking read(2) from C-code"
mktempdir() do d
    cfile=joinpath(d, "tmp.c")
    write(cfile, """
        #include <stdio.h>
        #include <stdlib.h>
        #include <unistd.h>
        #include <string.h>
        #include <sys/socket.h>

        #define LEN 100

        int main()
        {
            for(;;) {
                char buf[LEN];
                int n = read(0, buf, LEN);
                if (n == -1) {
                    perror("Exiting");
                    exit(0);
                }
                printf("%d:", n);
                for (int i = 0 ; i < n ; i++) {
                    if (buf[i] == '\\n') {
                        printf("\\\\n");
                    } else {
                        printf("%c", buf[i]);
                    }
                }
                printf("\\n");
                fflush(stdout);
            }
        }
        """)
    binfile=joinpath(d, "tmp")
    UnixIO.system("gcc -o $binfile $cfile")

    @info "readline() from pseudoterminal with fragmented writes."
    UnixIO.ptopen(`$binfile`; opts...) do cin, cout
        write(cin, "Hello1\n")
        @test readline(cout; keep=true) == "7:Hello1\\n\n"
        write(cin, "Hello2\n")
        @test readline(cout; keep=true) == "7:Hello2\\n\n"

        write(cin, "Hel")
        UnixIO.@cerr C.tcdrain(cin.fd)

        # Even after draining and waiting, the client read() has not returned
        # because of canonical mode.
        @test readline(cout; keep=true, timeout=1) == ""
        write(cin, "lo3\n")
        @test readline(cout; keep=true) == "7:Hello3\\n\n"

        write(cin, "Hello4\nHel")
        @test readline(cout; keep=true) == "7:Hello4\\n\n"

        @test readline(cout; keep=true, timeout=1) == ""
        write(cin, "lo5\n")
        @test readline(cout; keep=true) == "7:Hello5\\n\n"
    end

    @info "readline() from socket with fragmented writes."
    UnixIO.open(`$binfile`; opts...) do cin, cout
        write(cin, "Hello1\n")
        @test readline(cout; keep=true) == "7:Hello1\\n\n"
        write(cin, "Hello2\n")
        @test readline(cout; keep=true) == "7:Hello2\\n\n"
        write(cin, "Hel")
        @test readline(cout; keep=true, timeout=1) == "3:Hel\n"
        write(cin, "lo3\n")
        @test readline(cout; keep=true) == "4:lo3\\n\n"

        write(cin, "Hello4\nHel")
        @test readline(cout; keep=true) == "10:Hello4\\nHel\n"

        @test readline(cout; keep=true, timeout=1) == ""
        write(cin, "lo5\n")
        @test readline(cout; keep=true) == "4:lo5\\n\n"
    end
end

@info "readline() from pseudoterminal with fragmented writes."
mktempdir() do d
    cfile=joinpath(d, "tmp.c")
    write(cfile, """
        #include <stdlib.h>
        #include <unistd.h>
        #include <string.h>

        void write_s(char* s) {
            write(1, s, strlen(s));
        }

        int main()
        {
            sleep(1);
            write_s("Hel");
            sleep(1);
            write_s("lo1\\r\\n");
            write_s("Hello2\\r");
            sleep(1);
            write_s("\\nHello3");
            exit(0);
        }
        """)
    binfile=joinpath(d, "tmp")
    UnixIO.system("gcc -o $binfile $cfile")
    UnixIO.ptopen(`$binfile`; opts...) do cin, cout
        @test readline(cout; keep=true) == "Hello1\r\n"
        @test readline(cout; keep=true) == "Hello2\r\n"
        @test readline(cout; keep=true) == "Hello3"
    end
end

@info "readline() from pseudoterminal with fragmented writes (via bash)."
bash = "echo -n 'Hel'               ;" *
       "sleep 1                     ;" *
       "echo -n 'lo1\r\nHello2\r'   ;" *
       "sleep 1                     ;" *
       "echo -n '\nHello3'          ;"
UnixIO.ptopen(`bash -c "$bash"`; opts...) do cin, cout
    @test readline(cout) == "Hello1"
    @test readline(cout) == "Hello2"
    @test readline(cout) == "Hello3"
end

@info "readline() from socket with fragmented packets."
UnixIO.open(`cat`; opts...) do cin, cout

    write(cin, "Hello1\nHel")
    @test readline(cout) == "Hello1"

    # In socket mode, the partial line is in the fd buffer.
    @test String(take!(copy(cout.buffer))) == "Hel"

    write(cin, "lo2")
    r = Ref{Any}(nothing)
    @sync begin
        @async r[] = readline(cout)
        sleep(1)

        # readline() is waiting for the end of line.
        @test r[] == nothing
        close(cin)
    end

    # After close(), readline() compeltes.
    @test r[] == "Hello2"
end

@info "readline() from pseudoterminal with fragmented writes (via cat)."
UnixIO.ptopen(`cat`; opts...) do cin, cout
    @sync begin
        @async begin
            write(cin, "Hello0\nHello1\nHel")
            sleep(0.1)
            write(cin, "lo2\n")
            sleep(0.1)
            write(cin, "Hello3\nHel")
            write(cin, UInt8(C.CEOF)) # Terminate partial line.
            write(cin, UInt8(C.CEOF)) # Singal end of file
        end
        @test readline(cout; keep=true) == "Hello0\n"
        # In canonical mode, the partial "Hel" line is buffered by the OS
        # in cat's output buffer so only "Hello2\n" line is in the fd buffer.
        @test bytesavailable(cout.buffer) == 7
        @test String(take!(copy(cout.buffer))) == "Hello1\n"
        @test readline(cout; keep=true) == "Hello1\n"
        @test bytesavailable(cout.buffer) == 0
        @test readline(cout; keep=true) == "Hello2\n"
        @test readline(cout; keep=true) == "Hello3\n"
        # At end of transmission the partial line is returned:
        @test readline(cout; keep=true) == "Hel"
    end
end

@info "readline() from pseudoterminal non-canonlical mode."
UnixIO.ptopen(`cat`; opts...) do cin, cout

    UnixIO.tcsetattr(a -> a.c_lflag = 0, cout) # not C.ICANON

    @sync begin
        @async begin
            write(cin, "Hello0\nHello1\nHel")
            sleep(0.1)
            write(cin, "lo2")
            sleep(0.1)
            write(cin, UInt8(C.CEOF))
            sleep(0.1)
            write(cin, "Hello3\nHel")
            close(cin)
        end
        @test readline(cout; keep=true) == "Hello0\n"
        @test readline(cout; keep=true) == "Hello1\n"
        # In non canonical mode, the partial line is in the fd buffer:
        @test bytesavailable(cout.buffer) == 3
        # `CEOF` and sends a line with no "\n" (e.g. like the "bash$ " prompt)
        @test readline(cout; keep=true) == "Hello2\x04Hello3\n"
        # At end of transmission the partial line is returned:
        @test readline(cout; keep=true, timeout=1) == "Hel"
    end
end

hexdump_output = "00000000  48 65 6c 6c 6f 0a 48 65  6c 6c 6f 0a              |Hello.Hello.|\n0000000c\n"

@info "Short hexdump via socketpair."
p, out = UnixIO.open(`hexdump -C`; opts...) do cin, cout
    write(cin, "Hello\nHello\n")
    close(cin)
    read(cout, String)
end
@test out == hexdump_output

@info "Short hexdump via pseudoterminal."
p, out = UnixIO.ptopen(`hexdump -C`; opts...) do cin, cout
    println(cin, "Hello")
    println(cin, "Hello")
    close(cin)
    read(cout, String)
end
@test out == hexdump_output

@info "Short pseudoterminal hexdump, read before writet."
p, out = UnixIO.ptopen(`hexdump -C`; opts...) do cin, cout
    @sync begin
        @async begin
            println(cin, "Hello")
            println(cin, "Hello")
            close(cin)
        end
        read(cout, String)
    end
end
@test out == hexdump_output

@info "Short pseudoterminal hexdump, interleaved read/write."
p, out = UnixIO.ptopen(`hexdump -C`; opts...) do cin, cout
    @sync begin
        println(cin, "Hello")
        @async begin
            println(cin, "Hello")
            close(cin)
        end
        read(cout, String)
    end
end
@test out == hexdump_output

@info "Interactive bash script via pseudoterminal."
#env = merge(ENV, Dict("TERM" => "dumb"))
UnixIO.ptopen(`bash`; opts...) do cin, cout
    ptdump(`bash`, cin, cout) do fin, fout
        # Wait for "bash-3.2$" prompt.
        while !contains(fin(), r"\$") end

        fout(raw"""i=1""" * "\n")                 ;
        @test contains(fin(), r"\$")
        fout(raw"""while [ $i -lt 11 ]""" * "\n") ; @test fin() == "> "
        fout(raw"""do""" * "\n")                  ; @test fin() == "> "
        fout(raw"""    sleep 0.2""" * "\n")       ; @test fin() == "> "
        fout(raw"""    echo "COUNT$i" """ * "\n") ; @test fin() == "> "
        fout(raw"""    i=$(($i+1))""" * "\n")     ; @test fin() == "> "
        fout(raw"""done""" * "\n")

        for i in 1:10
            @test fin() == "COUNT$i\n"
        end
        @test contains(fin(), r"\$")
    end
end

@info "Interactive julia script via pseudoterminal."
env = merge(ENV, Dict("TERM" => "dumb"))

UnixIO.ptopen(`$(Base.julia_cmd()) --color=yes`; env=env) do cin, cout

    ptdump(Base.julia_cmd(), cin, cout) do fin, fout

        # Wait for "julia>" prompt.
        while !contains(fin(), r"julia> ") end

        fout(raw"""ENV["TERM"]""" * "\n")
        @test fin() == "\"dumb\"\n"
        @test fin() == "\n"
        @test fin() == "julia> "

        fout(raw"""i=1""" * "\n")
        @test fin() == "1\n"
        @test fin() == "\n"
        @test fin() == "julia> "

        fout(raw"""while i < 11""" * "\n")
        fout(raw"""    sleep(0.2)""" * "\n")
        fout(raw"""    println("COUNT$i")""" * "\n")
        fout(raw"""    i += 1""" * "\n")
        fout(raw"""end""" * "\n")

        for i in 1:10
            @test fin() == "COUNT$i\n"
        end
        @test fin() == "\n"
        @test fin() == "julia> "
    end
end

end # @testset "Pseudoterminal Tests"


text:
@testset "LTI system kernels tests" begin
    N = 10
    types = [Int32, Int64, Float32, Float64, ComplexF32, ComplexF64]
    atol = 1e-2

    @testset "differentiator tests" begin
        Δt = 0.1
        s = Romeo.LTI.Diff{Float64}(1)

        k = Romeo.LTI.kernel(s, N, Δt)
        @test k[1] == 1/Δt || "The first element of the differentiation kernel should be equal to 1/Δt"
        @test k[2] == -1/Δt || "The second element of the differentiation kernel should be -1/Δt"
        for i = 3:N
            @test k[i] == 0 || "All elements with index ≥3 of the differentiation kernel should be zero."
        end # for

        ys = Romeo.LTI.step_resp(s, N, Δt)
        @test ys[1] == 1/Δt || "The first element of the step response of a differentiator should be 1/Δt."
        for i = 2:N
            @test ys[i] == 0 || "All elements of the step response of a differentiator should be zero, except the first one."
        end # for
    end # testset

    @testset "identity operator tests" begin
        Δt = 0.1
        s = Romeo.LTI.Diff{Float64}(0)

        k = Romeo.LTI.kernel(s, N, Δt)
        @test k[1] == 1 || "The first element of the identity kernel should be equal to 1"
        for i = 2:N
            @test k[i] == 0 || "All elements with index ≥2 of the identity kernel should be zero."
        end # for

        ys = Romeo.LTI.step_resp(s, N, Δt)
        for i = 1:N
            @test ys[i] == 1 || "All elements of the step response of the identity kernel should be equal to 1."
        end # for
    end # testset

    @testset "integrator tests" begin
        Δt = 1
        sinv = Romeo.LTI.Diff{Float64}(-1)
        k = Romeo.LTI.kernel(sinv, N, Δt)
        for i = 1:N
            @test k[i] == Δt || "All elements of the integration kernel should be equal to Δt."
        end # for

        ys = Romeo.LTI.step_resp(sinv, N, Δt)
        for i = 1:N
            @test ys[i] == i*Δt || "All elements of the integration kernel should be equal to i*Δt."
        end # for
    end # testset

    @testset "fractional integrator tests" begin
        Δt = 1
        int_0_5 = Romeo.LTI.Diff{Float64}(-0.5)
        int_1_5 = Romeo.LTI.Diff{Float64}(-0.5)

        ys = Romeo.LTI.impulse_resp(int_0_5, N, Δt)
        yi = Romeo.LTI.impulse_resp(int_1_5, N, Δt)

        @test all(ys .== yi) || "Step response of half integrator should match impulse response of 1.5 integrator."
    end # testset

end # testset

ext: jl


text:
module TestApproxRidge
    using StreamStats
    using Distributions
    using Base.Test

    for n in rand(100_000:1_000_000, 10)
        p = rand(1:100)
        β₀ = randn()
        β = randn(p)
        xs = randn(p, n)
        ys = Float64[β₀ + dot(β, xs[:, i]) + randn() for i in 1:n]
        stat = StreamStats.ApproxRidge(p, 1.0)
        for i in 1:n
            update!(stat, xs[:, i], ys[i])
        end
        @test cor(vcat(β₀, β), state(stat)) > 0.99
        # TODO: Add real tests
        # l2_err = norm(vcat(β₀, β) - state(stat))
        # @test l2_err < 0.05 * p
        # @test nobs(stat) == n
    end
end


text:
module weird_timing_examples

using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforLinear: FEMMDeforLinear, lumpedmass
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.MatDeforNeohookeanNaiveModule: MatDeforNeohookeanNaive
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness, geostiffness, nzebcloads, restoringforce, estimatestablestep
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearExplModule: FEMMDeforNonlinearExpl
using FinEtoolsDeforNonlinear.AssemblyModule: SysvecAssemblerOpt
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots

mr = DeforModelRed3D
E, nu = 7.0*phun("MPa"), 0.3
mass_density = 1000.0*phun("kg/m^3")
L= 6/2*phun("mm");
H = 2/2*phun("mm");
W = 2/2*phun("mm");
tmag = 0.1*phun("MPa");# Magnitude of the traction
tolerance = W / 1000
traction_vector = [0.0, 0.0, -tmag]
# A million-element mesh
nL, nW, nH = 160, 80, 80
tend = 0.00075e-3
# Much smaller mesh
nL, nW, nH = 80, 40, 40
tend = 0.001e-3

function neohookean_h8()
	timing = time()

	m = MatDeforNeohookean(mr, mass_density, E, nu)
	fens, fes = H8block(L, W, H, nL, nW, nH)
    @info count(fens), count(fes)
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3))

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    for i in 1:3
        setebc!(u, l1, true, i, 0.0)
    end
    applyebc!(u)
    numberdofs!(u)

    bfes = meshboundary(fes)
    el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    function setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
        v .= 1.0 .* traction_vector
        return v
    end
    fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
    eL1femm =  FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2)))
    FL = distribloads(eL1femm, geom, u, fi, 2);

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

    femm = FEMMDeforNonlinearExpl(mr, IntegDomain(fes, GaussRule(3, 2)), m)
    femm = associategeometry!(femm, geom)

    Ux = FFlt[]; ts = FFlt[]
    function increment_observer(step, t, un1)
        if step == 1 || rem(step, 10) == 0
            println("$(step) steps")
            push!(Ux, mean(un1.values[movel1,3]));
            push!(ts, t);
        end
    end

    un1 = deepcopy(u)
    un = deepcopy(u)
    vn1 = deepcopy(u)

    stabldt = estimatestablestep(femm, geom);
    M = lumpedmass(femm, geom, un1)
    invMv = [1.0 / M[idx, idx] for idx in 1:size(M, 1)]

    dtn = 0.8 * stabldt
    tn = 0.0
    # Initial displacement, velocity, and acceleration.
    Un = gathersysvec(un1);
    Vn = gathersysvec(vn1);
    # The acceleration will be computed from the initial loads.
    An = deepcopy(Vn);
    Fn = deepcopy(Vn);
    # Temporary vectors
    Un1 = deepcopy(Un);
    Vn1 = deepcopy(Vn);
    An1 = deepcopy(An);
    # Global vector assembler
    assembler = SysvecAssemblerOpt(deepcopy(Un1), length(Un1))
    step = 0;
    while tn < tend
        step = step + 1      # Step counter
        fill!(Fn, 0.0) # Zero out the load
        Fn .= Fn .+ FL # Add on the time-independent load vector
        # If this is the first step compute the initial acceleration.
        if step == 1
            An .= invMv .* Fn;
            # An .= M \ Fn;
        end
        # Update the displacements
        @. Un1 = Un + dtn*Vn + ((dtn^2)/2)*An;# displacement update
        scattersysvec!(un1, Un1);
        # Add the restoring forces, starting from the time-independent load.
        tim = time()
        Fn .+= restoringforce(femm, assembler, geom, un1, un, tn, dtn, true)
        println("$(time() - tim)")
        # Compute the new acceleration.
        An1 .= invMv .* Fn;
        # An1 .= M \ Fn;
        # Update the velocity
        @. Vn1 = Vn + (dtn/2) * (An+An1);
        # Bring the the displacement and velocity fields up to date
        # scattersysvec!(vn1, Vn1);
        # Switch the temporaries for the next step.
        (Un, Un1) = (Un1, Un);
        (Vn, Vn1) = (Vn1, Vn);
        (An, An1) = (An1, An);
        (un, un1) = (un1, un);
        if (tn == tend)   # are we at the end?
            break;
        end
        if (tn+dtn > tend) # If the next step would take us beyond the end,
            dtn = tend-tn; # hit the end precisely.
        end
        tn = tn+dtn;
        increment_observer(step, tn, un1);
    end

	timing = time() - timing

    @show ts, Ux / phun("mm")
    # pl = lineplot(ts, Ux / phun("mm"), canvas = DotCanvas)
    # display(pl)

    # vtkexportmesh("neohookean_h8.vtk", fens, fes; vectors = [("u", un1.values)])

    println("$step steps done in $(timing) seconds")
    mus = timing / count(fes) / step * 1.0e6
    println("$(round(100 * mus) / 100) microseconds per element in one time step")

    true
end # function neohookean_h8

# For threaded execution we will need a per-thread buffer:
# Each thread has its own machine, its own assembler, and its own assembled
# vector. This is to avoid memory contention.
struct ThreadBuffer
	femm # finite element machine, with private copy of finite elements and the material
	assembler # assembler: each thread needs to have its own in order not to write into one the same
end

function neohookean_h8_thr(NTHREADS)
	timing = time()

    m = MatDeforNeohookean(mr, mass_density, E, nu)
    fens, fes = H8block(L, W, H, nL, nW, nH)
    @info "Mesh of $(count(fens)) nodes, $(count(fes)) elements"
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3))

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    for i in 1:3
        setebc!(u, l1, true, i, 0.0)
    end
    applyebc!(u)
    numberdofs!(u)

    bfes = meshboundary(fes)
    el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    function setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
        v .= 1.0 .* traction_vector
        return v
    end
    fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
    eL1femm =  FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2)))
    FL = distribloads(eL1femm, geom, u, fi, 2);

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

    # Now we prepare  the assembly for threaded execution.
    @show nth = NTHREADS #Base.Threads.nthreads()
    chunk = Int(floor(count(fes) / nth))
    threadbuffs = ThreadBuffer[];
    for th in 1:nth
    	# This thread will work on a subset of the mesh.
    	meshrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:count(fes))
    	feschnk = deepcopy(subset(fes, collect(meshrange)))
    	material = deepcopy(m)
    	femm = FEMMDeforNonlinearExpl(mr, IntegDomain(feschnk, GaussRule(3, 2)), material)
    	femm = associategeometry!(femm, geom)
    	# It will use its own assembler in order to avoid memory contention.
    	assembler = SysvecAssemblerOpt(fill(0.0, u.nfreedofs), u.nfreedofs)
    	push!(threadbuffs, ThreadBuffer(femm, assembler));
    end

    Ux = FFlt[]; ts = FFlt[]
    function increment_observer(step, t, un1)
        if step == 1 || rem(step, 3) == 0
            println("$(step) steps")
            push!(Ux, mean(un1.values[movel1,3]));
            push!(ts, t);
        end
    end

    # Fields for the current, last displacements
    un1 = deepcopy(u)
    un = deepcopy(u)

    # Assemble the mass matrix for the entire mesh
    femm = FEMMDeforNonlinearExpl(mr, IntegDomain(fes, GaussRule(3, 2)), m)
    stabldt = estimatestablestep(femm, geom);
    M = lumpedmass(femm, geom, un1)
    invMv = [1.0 / M[idx, idx] for idx in 1:size(M, 1)]
    femm = nothing

    # Prepare the time stepping
    dtn = 0.8 * stabldt
    tn = 0.0
    # Initial displacement, velocity, and acceleration.
    Un = gathersysvec(un1);
    Vn = gathersysvec(un1);
    # The acceleration will be computed from the initial loads.
    An = deepcopy(Vn);
    Fn = deepcopy(Vn);
    # Temporary vectors
    Un1 = deepcopy(Un);
    Vn1 = deepcopy(Vn);
    An1 = deepcopy(An);
    timfo = []
    timco = []
    timto = []
    step = 0;
    while tn < tend
        step = step + 1      # Step counter
        fill!(Fn, 0.0) # Zero out the load
        Fn .= Fn .+ FL # Add on the time-independent load vector
        # If this is the first step compute the initial acceleration.
        if step == 1
            An .= invMv .* Fn;
        end
        # Update the displacements
        @. Un1 = Un + dtn*Vn + ((dtn^2)/2)*An;# displacement update
        scattersysvec!(un1, Un1);
        # Add the restoring forces, starting from the time-independent load.
        tim = time()
        tim1  = time()
        tasks = [];
        for th in 1:length(threadbuffs)
        	push!(tasks, Threads.@spawn begin
        		fill!(threadbuffs[th].assembler.F_buffer, 0.0)
        		# now add the restoring force from the subset of the mesh handled by this thread
        		restoringforce(threadbuffs[th].femm, threadbuffs[th].assembler, geom, un1, un, tn, dtn, true)
        	end);
        end
        push!(timfo, time() - tim1)
        # Wait for the threads to finish, and then add the force from the thread to the global force vector
        tim1  = time()
        for th in 1:length(tasks)
        	Threads.wait(tasks[th]);
        	Fn .+= threadbuffs[th].assembler.F_buffer
        end
        push!(timco, time() - tim1)
        push!(timto, time() - tim)
        # Compute the new acceleration.
        An1 .= invMv .* Fn;
        # An1 .= M \ Fn;
        # Update the velocity
        @. Vn1 = Vn + (dtn/2) * (An+An1);
        # Bring the the displacement and velocity fields up to date
        # scattersysvec!(vn1, Vn1);
        # Switch the temporaries for the next step.
        (Un, Un1) = (Un1, Un);
        (Vn, Vn1) = (Vn1, Vn);
        (An, An1) = (An1, An);
        (un, un1) = (un1, un);
        if (tn == tend)   # are we at the end?
            break;
        end
        if (tn+dtn > tend) # If the next step would take us beyond the end,
            dtn = tend-tn; # hit the end precisely.
        end
        tn = tn+dtn;
        increment_observer(step, tn, un1);
    end

	timing = time() - timing

    ts, Ux / phun("mm")
    @show mean(timto), nth
    # pl = lineplot(ts, Ux / phun("mm"), canvas = DotCanvas)
    # display(pl)

    # vtkexportmesh("neohookean_h8_thr.vtk", fens, fes; vectors = [("u", un1.values)])

    println("$step steps done in $(timing) seconds")
    mus = timing / count(fes) / step * 1.0e6
    println("$(round(100 * mus) / 100) microseconds per element in one time step")

    true
end # function neohookean_h8_thr

function neohookean_h8_thr_2(NTHREADS)
	timing = time()

    m = MatDeforNeohookean(mr, mass_density, E, nu)
    	fens, fes = H8block(L, W, H, nL, nW, nH)
    @info "Mesh of $(count(fens)) nodes, $(count(fes)) elements"
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3))

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    for i in 1:3
        setebc!(u, l1, true, i, 0.0)
    end
    applyebc!(u)
    numberdofs!(u)

    bfes = meshboundary(fes)
    el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    function setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
        v .= 1.0 .* traction_vector
        return v
    end
    fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
    eL1femm =  FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2)))
    FL = distribloads(eL1femm, geom, u, fi, 2);

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

    # Now we prepare  the assembly for threaded execution.
    @show nth = NTHREADS #Base.Threads.nthreads()
    chunk = Int(floor(count(fes) / nth))
    threadbuffs = ThreadBuffer[];
    for th in 1:nth
    	# This thread will work on a subset of the mesh.
    	meshrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:count(fes))
    	feschnk = deepcopy(subset(fes, collect(meshrange)))
    	material = deepcopy(m)
    	femm = FEMMDeforNonlinearExpl(mr, IntegDomain(feschnk, GaussRule(3, 2)), material)
    	femm = associategeometry!(femm, geom)
    	# It will use its own assembler in order to avoid memory contention.
    	assembler = SysvecAssemblerOpt(fill(0.0, u.nfreedofs), u.nfreedofs)
    	push!(threadbuffs, ThreadBuffer(femm, assembler));
    end

    Ux = FFlt[]; ts = FFlt[]
    function increment_observer(step, t, un1)
    	if step == 1 || rem(step, 3) == 0
    		println("$(step) steps")
    		push!(Ux, mean(un1.values[movel1,3]));
    		push!(ts, t);
    	end
    end

    # Fields for the current, last displacements
    un1 = deepcopy(u)
    un = deepcopy(u)

    # Assemble the mass matrix for the entire mesh
    femm = FEMMDeforNonlinearExpl(mr, IntegDomain(fes, GaussRule(3, 2)), m)
    stabldt = estimatestablestep(femm, geom);
    M = lumpedmass(femm, geom, un1)
    invMv = [1.0 / M[idx, idx] for idx in 1:size(M, 1)]
    femm = nothing

    # Prepare the time stepping
    dtn = 0.8 * stabldt
    tn = 0.0
    # Initial displacement, velocity, and acceleration.
    Un = gathersysvec(un1);
    Vn = gathersysvec(un1);
    # The acceleration will be computed from the initial loads.
    An = deepcopy(Vn);
    Fn = deepcopy(Vn);
    # Temporary vectors
    Un1 = deepcopy(Un);
    Vn1 = deepcopy(Vn);
    An1 = deepcopy(An);
    timfo = []
    timco = []
    timto = []
    step = 0;
    while tn < tend
        step = step + 1      # Step counter
        fill!(Fn, 0.0) # Zero out the load
        Fn .= Fn .+ FL # Add on the time-independent load vector
        # If this is the first step compute the initial acceleration.
        if step == 1
            An .= invMv .* Fn;
        end
        # Update the displacements
        @. Un1 = Un + dtn*Vn + ((dtn^2)/2)*An;# displacement update
        scattersysvec!(un1, Un1);
        # Add the restoring forces, starting from the time-independent load.
        tim = time()
        tim1  = time()
        Threads.@threads for th in 1:length(threadbuffs)
        begin
	        fill!(threadbuffs[th].assembler.F_buffer, 0.0)
	        # now add the restoring force from the subset of the mesh handled by this thread
	        restoringforce(threadbuffs[th].femm, threadbuffs[th].assembler, geom, un1, un, tn, dtn, true)
	    end
	    end
        push!(timfo, time() - tim1)
        # Wait for the threads to finish, and then add the force from the thread to the global force vector
        tim1  = time()
        for th in 1:length(threadbuffs)
        	Fn .+= threadbuffs[th].assembler.F_buffer
        end
        push!(timco, time() - tim1)
        push!(timto, time() - tim)
        # Compute the new acceleration.
        An1 .= invMv .* Fn;
        # An1 .= M \ Fn;
        # Update the velocity
        @. Vn1 = Vn + (dtn/2) * (An+An1);
        # Bring the the displacement and velocity fields up to date
        # scattersysvec!(vn1, Vn1);
        # Switch the temporaries for the next step.
        (Un, Un1) = (Un1, Un);
        (Vn, Vn1) = (Vn1, Vn);
        (An, An1) = (An1, An);
        (un, un1) = (un1, un);
        if (tn == tend)   # are we at the end?
            break;
        end
        if (tn+dtn > tend) # If the next step would take us beyond the end,
            dtn = tend-tn; # hit the end precisely.
        end
        tn = tn+dtn;
        increment_observer(step, tn, un1);
    end

	timing = time() - timing

    @show ts, Ux / phun("mm")
    @show mean(timto), nth
    # pl = lineplot(ts, Ux / phun("mm"), canvas = DotCanvas)
    # display(pl)

    # vtkexportmesh("neohookean_h8_thr.vtk", fens, fes; vectors = [("u", un1.values)])

    println("$step steps done in $(timing) seconds")
    mus = timing / count(fes) / step * 1.0e6
    println("$(round(100 * mus) / 100) microseconds per element in one time step")

    true
end # function neohookean_h8_thr

function allrun()
    println("#####################################################")
    println("# neohookean_h8 ")
    neohookean_h8()             # weird_timing_examples.neohookean_h8()
    println("#####################################################")
    println("# neohookean_h8_thr ")
    neohookean_h8_thr()             # weird_timing_examples.neohookean_h8_thr()
    println("#####################################################")
    println("# neohookean_h8_thr_2 ")
    neohookean_h8_thr_2()
    return true
end # function allrun

end # module weird_timing_examples


text:
using Documenter, Slacker

makedocs(;
    modules=[Slacker],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
    ],
    repo="https://github.com/Sumidu/slacker.jl/blob/{commit}{path}#L{line}",
    sitename="slacker.jl",
    authors="André Calero Valdez",
    assets=String[],
)

deploydocs(;
    repo="github.com/Sumidu/slacker.jl",
)


text:
module SSLTestOverloads
using SimpleSDMLayers
using Test

M = rand(Bool, (5, 10))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)

for i in eachindex(M)
   @test S[i] == M[i]
end
for i in eachindex(S)
   @test S[i] == M[i]
end

@test typeof(S[0.2, 0.6]) == eltype(M)
@test isnothing(S[1.2, 0.3])
@test isnothing(S[1.2, 1.3])
@test isnothing(S[0.2, 1.3])

@test typeof(S[1:2, 5:7]) == typeof(S)

@test typeof(S[left=0.2, right=0.6, bottom=0.5, top=1.0]) == typeof(S)
@test S[left=0.2, right=0.6, bottom=0.5, top=1.0].top ≈ 1.0
@test S[left=0.2, right=0.6, bottom=0.5, top=1.0].bottom ≈ 0.4
@test S[left=0.2, right=0.6, bottom=0.5, top=1.0].right ≈ 0.6
@test S[left=0.19, right=0.6, bottom=0.5, top=1.0].left <= 0.2

@test typeof(S[left=0.2, bottom=0.5]) == typeof(S)
@test S[left=0.2, bottom=0.5].top ≈ 1.0
@test S[left=0.2, bottom=0.5].bottom ≈ 0.4
@test S[left=0.2, bottom=0.5].right ≈ 1.0
@test S[left=0.2, bottom=0.5].left ≈ 0.2

C = (left=0.2, bottom=0.5)
@test typeof(S[C]) == typeof(S)

Y = SimpleSDMResponse(zeros(Float64, (5,5)), 0.0, 1.0, 0.0, 1.0)
Y[0.1,0.1] = 0.2
@test Y[0.1,0.1] == 0.2

Z = convert(SimpleSDMPredictor, Y)
Y[0.1,0.1] = 4.0
@test Z[0.1,0.1] != Y[0.1,0.1]
@test Z[0.1,0.1] == 0.2

Z = convert(SimpleSDMPredictor, Y)
V = collect(Z)
@test typeof(V) == Vector{Float64}

# hcat / vcat
l1 = worldclim(1, left=0.0, right=10.0, bottom=0.0, top=10.0)
l2 = worldclim(1, left=0.0, right=10.0, bottom=10.0, top=20.0)
l3 = worldclim(1, left=10.0, right=20.0, bottom=0.0, top=10.0)
l4 = worldclim(1, left=10.0, right=20.0, bottom=10.0, top=20.0)

ml1 = hcat(l1, l3)
vl1 = vcat(l1, l2)
ml2 = hcat(l2, l4)
vl2 = vcat(l3, l4)

@test all(vcat(ml1, ml2).grid == hcat(vl1, vl2).grid)

# typed similar
c2 = similar(Bool, l1)
@test eltype(c2) == Bool

end


text:
using JLD, FileIO
using TikzGraphs
using JSON
using ProgressMeter
include("./network_library.jl")

function make_figs(experiment_dict::Dict,f_loc::String,i_loc::String)
    for expr in keys(experiment_dict[:conditions])
        net_type = experiment_dict[:name]
        fn = experiment_dict[:conditions][expr][:fname]
        fldr = experiment_dict[:conditions][expr][:fldr]

        println("Processing: $fn")
        println("running condition '$expr' data")

        f = joinpath(fldr,fn)

        info("# Running $(f)")
        df = JLD.load("$f.jld")
        d = df["problem_dict"]
        g = df["global_rwd_range"]

        info("Processing Networks")
        p = Progress(length(d),barglyphs=BarGlyphs("[=> ]"),color=:white)
        j = 0
        for net in d
            net_num = net[1]
            net = net[2]
            e_start = net[:evader_start]
            p_start = net[:pursuer_start]
            net_fname = "$(net_type)_$(expr)_net_$net_num"
            f_type = :svg
            display_network(net[:graph],evader_locs=[e_start],pursuer_locs=[p_start],scale=1.0,
                            fname=joinpath(i_loc,net_fname),ftype=f_type)

            net[:img_fldr] = i_loc
            net[:img_name] = net_fname * "." * String(f_type)
            j += 1
            ProgressMeter.next!(p; showvalues = [(:network,"$j of $(length(d))")])
        end

        # re-write original file with "r_dist" object added to each network
        info("Writing to $f")
        jldopen("$f.jld", "w") do file
            write(file, "problem_dict", d)
            write(file, "global_rwd_range", g)
        end
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using Pkg

# Pkg.add("BinaryBuilder")
Pkg.add("IJulia")
Pkg.add("NamedArrays")
Pkg.add("Feather")
Pkg.add("DataFrames")
Pkg.add("Conda")
# Problems with Plots > GR > GR_jll > Qt5Base:
# - dylib not found
# - if renamed 'binary' to .dylib >> issue with _NSAppearanceNameDarkAqua not found (this is only avail. on 10.14!))
# cf. here: https://discourse.julialang.org/t/problem-with-qt-dependency-of-plots-jl/58765
# cf. here: https://github.com/jheinen/GR.jl/issues/392
# v1.1.2 is causing the dylib not found error (and subsequently the _NSAppearanceNameDarkAqua when renaming the framework binaries to dylib)
# v1.1.1 is causing the _NSAppearanceNameDarkAqua not found error
ENV["JULIA_GR_PROVIDER"] = "GR"
Pkg.add(Pkg.PackageSpec(;name="Plots", version="1.1.0"))

Pkg.add("WebIO")
Pkg.add("Interact")

using Plots
using IJulia
using WebIO

# precompile the packages ?
Pkg.precompile()

max_stars_repo_name: benbenz/jupyterlab-econ-desktop

max_issues_repo_name: benbenz/jupyterlab-econ-desktop

max_forks_repo_name: benbenz/jupyterlab-econ-desktop


text:
"""
    PowerNetwork <: AbstractPowerSystem

This structure contains constant parameters that define the topology and
physics of the power network.

The object `PowerNetwork` uses its own contiguous indexing for the buses.
The indexing is independent from those specified in the Matpower or the
PSSE input file. However, a correspondence between the two indexing
(Input indexing *to* `PowerNetwork` indexing) is stored inside the
attribute `bus_to_indexes`.

## Note
The object `PowerNetwork` is created in the main memory.
Use a `AbstractFormulation` object to move data to the target device.

"""
struct PowerNetwork <: AbstractPowerSystem
    vbus::Vector{Complex{Float64}}
    Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}
    data::Dict{String,Array}

    nbus::Int64
    ngen::Int64

    bustype::Vector{Int64}
    bus_to_indexes::Dict{Int, Int}
    ref::Vector{Int64}
    pv::Vector{Int64}
    pq::Vector{Int64}

    sbus::Vector{Complex{Float64}}
    sload::Vector{Complex{Float64}}

    function PowerNetwork(datafile::String, data_format::Int64)

        if data_format == 0
            println("Reading PSSE format")
            data_raw = ParsePSSE.parse_raw(datafile)
            data, bus_id_to_indexes = ParsePSSE.raw_to_exapf(data_raw)
        elseif data_format == 1
            data_mat = ParseMAT.parse_mat(datafile)
            data, bus_id_to_indexes = ParseMAT.mat_to_exapf(data_mat)
        end
        # Parsed data indexes
        BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
        LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()

        # retrive required data
        bus = data["bus"]
        gen = data["gen"]
        SBASE = data["baseMVA"][1]

        # size of the system
        nbus = size(bus, 1)
        ngen = size(gen, 1)

        # obtain V0 from raw data
        vbus = zeros(Complex{Float64}, nbus)
        for i in 1:nbus
            vbus[i] = bus[i, VM]*exp(1im * pi/180 * bus[i, VA])
        end

        # form Y matrix
        Ybus = makeYbus(data, bus_id_to_indexes)

        # bus type indexing
        ref, pv, pq, bustype = bustypeindex(bus, gen, bus_id_to_indexes)

        sbus, sload = assembleSbus(gen, bus, SBASE, bus_id_to_indexes)

        new(vbus, Ybus, data, nbus, ngen, bustype, bus_id_to_indexes, ref, pv, pq, sbus, sload)
    end
end
function PowerNetwork(datafile::String)
    if endswith(datafile, ".raw")
        data_format = 0
    elseif endswith(datafile, ".m")
        data_format = 1
    else
        error("Unsupported format in file $(datafile): supported extensions are " *
              "Matpower (.m) or PSSE (.raw)")
    end
    return PowerNetwork(datafile, data_format)
end

# Getters
## Network attributes
get(pf::PowerNetwork, ::NumberOfBuses) = pf.nbus
get(pf::PowerNetwork, ::NumberOfLines) = size(pf.data["branch"], 2)
get(pf::PowerNetwork, ::NumberOfGenerators) = pf.ngen
get(pf::PowerNetwork, ::NumberOfPVBuses) = length(pf.pv)
get(pf::PowerNetwork, ::NumberOfPQBuses) = length(pf.pq)
get(pf::PowerNetwork, ::NumberOfSlackBuses) = length(pf.ref)

## Indexing
function get(pf::PowerNetwork, ::GeneratorIndexes)
    GEN_BUS = IndexSet.idx_gen()[1]
    gens = pf.data["gen"]
    ngens = size(gens, 1)
    # Create array on host memory
    indexing = zeros(Int, ngens)
    # Here, we keep the same ordering as specified in Matpower.
    for i in 1:ngens
        indexing[i] = pf.bus_to_indexes[gens[i, GEN_BUS]]
    end
    return indexing
end
get(pf::PowerNetwork, ::PVIndexes) = pf.pv
get(pf::PowerNetwork, ::PQIndexes) = pf.pq
get(pf::PowerNetwork, ::SlackIndexes) = pf.ref

# Pretty printing
function Base.show(io::IO, pf::PowerNetwork)
    println("Power Network characteristics:")
    @printf("\tBuses: %d. Slack: %d. PV: %d. PQ: %d\n", pf.nbus, length(pf.ref),
            length(pf.pv), length(pf.pq))
    println("\tGenerators: ", pf.ngen, ".")
    # Print system status
    @printf("\t==============================================\n")
    @printf("\tBUS \t TYPE \t VMAG \t VANG \t P \t Q\n")
    @printf("\t==============================================\n")

    for i=1:pf.nbus
        type = pf.bustype[i]
        vmag = abs(pf.vbus[i])
        vang = angle(pf.vbus[i])*(180.0/pi)
        pinj = real(pf.sbus[i])
        qinj = imag(pf.sbus[i])
        @printf("\t%d \t  %d \t %1.3f\t%3.2f\t%3.3f\t%3.3f\n", i,
                type, vmag, vang, pinj, qinj)
    end
end

# Some utils function

function bounds(pf::PowerNetwork, ::Buses, ::VoltageMagnitude)
    BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
    LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()

    bus = pf.data["bus"]
    v_min = convert.(Float64, bus[:, VMIN])
    v_max = convert.(Float64, bus[:, VMAX])
    return v_min, v_max
end

function bounds(pf::PowerNetwork, ::Generator, ::ActivePower)
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()

    gens = pf.data["gen"]
    baseMVA = pf.data["baseMVA"][1]

    p_min = convert.(Float64, gens[:, PMIN] / baseMVA)
    p_max = convert.(Float64, gens[:, PMAX] / baseMVA)
    return p_min, p_max
end

function bounds(pf::PowerNetwork, ::Generator, ::ReactivePower)
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()

    gens = pf.data["gen"]
    baseMVA = pf.data["baseMVA"][1]

    q_min = convert.(Float64, gens[:, QMIN] / baseMVA)
    q_max = convert.(Float64, gens[:, QMAX] / baseMVA)
    return q_min, q_max
end

"""
    get_costs_coefficients(pf::PowerNetwork)

Return coefficients for costs function.

"""
function get_costs_coefficients(pf::PowerNetwork)
    # indexes
    BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
    LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()
    MODEL, STARTUP, SHUTDOWN, NCOST, COST = IndexSet.idx_cost()

    ref = pf.ref
    pv = pf.pv
    pq = pf.pq
    b2i = pf.bus_to_indexes

    # Matpower assumes gens are ordered. Generator in row i has its cost on row i
    # of the cost table.
    gens = pf.data["gen"]
    baseMVA = pf.data["baseMVA"][1]
    bus = pf.data["bus"]
    ngens = size(gens)[1]
    nbus = size(bus)[1]

    # initialize cost
    # store coefficients in a Float64 array, with 4 columns:
    # - 1st column: bus type
    # - 2nd column: constant coefficient c0
    # - 3rd column: coefficient c1
    # - 4th column: coefficient c2
    coefficients = zeros(ngens, 4)
    # If cost is not specified, we return the array coefficients as is
    if !haskey(pf.data, "cost")
        @warn("PowerSystem: cost is not specified in PowerNetwork dataset")
        return coefficients
    end

    cost_data = pf.data["cost"]
    # iterate generators and check if pv or ref.
    for i = 1:ngens
        # only 2nd degree polynomial implemented for now.
        @assert cost_data[i, MODEL] == 2
        @assert cost_data[i, NCOST] == 3
        genbus = b2i[gens[i, GEN_BUS]]
        bustype = bus[genbus, BUS_TYPE]

        # polynomial coefficients
        # TODO: currently scale by baseMVA. Is it a good idea?
        c0 = cost_data[i, COST][3]
        c1 = cost_data[i, COST][2] * baseMVA
        c2 = cost_data[i, COST][1] * baseMVA^2
        coefficients[i, :] .= (bustype, c0, c1, c2)
    end
    return coefficients
end


text:
_length(rws) = length(collect(KnuthBendix.rules(rws)))

@testset "KBS1 examples" begin
    R = RWS_Example_5_1()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs1)
    @test _length(rws) == 8
    @test collect(KnuthBendix.rules(R))[1:5] == collect(KnuthBendix.rules(rws))[1:5]

    R = RWS_Example_5_2()
    rws = @test_logs (:warn,) KnuthBendix.knuthbendix(R, maxrules=100, implementation=:naive_kbs1)
    @test _length(rws) > 50 # there could be less rules that 100 in the irreducible rws

    R = RWS_Example_5_3()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs1)
    @test _length(rws) == 6
    a,b = Word.([i] for i in 1:2)
    ε = one(a)
    @test collect(KnuthBendix.rules(rws)) == KnuthBendix.Rule.([a^2=>ε, b^3=>ε,
        (b*a)^2=>a*b^2, (a*b)^2=>b^2*a, a*b^2*a=>b*a*b, b^2*a*b^2=>a*b*a])

    R = RWS_Example_5_4()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs1)
    @test _length(rws) == 11

    R = RWS_Example_5_5()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs1)
    @test _length(rws) == 18

    R = RWS_Example_6_4()
    rws = KnuthBendix.knuthbendix(R, maxrules=100, implementation=:naive_kbs1)
    @test _length(rws) == 40

    R = RWS_Example_6_5()
    rws = KnuthBendix.knuthbendix(R, maxrules=100, implementation=:naive_kbs1)
    @test _length(rws) == 40

    R = RWS_Closed_Orientable_Surface(3)
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs1)
    @test _length(rws) == 16
end

function test_kbs2_methods(R, methods, len)
    rwses = [KnuthBendix.knuthbendix(R, implementation=m) for m in methods]
    lengths = _length.(rwses)
    @test all(==(len), lengths)
end

@testset "KBS2 examples" begin
    R = RWS_Example_5_1()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs2)
    @test _length(rws) == 8
    @test Set(collect(KnuthBendix.rules(R))[1:5]) ==
        Set(collect(KnuthBendix.rules(rws))[1:5])

    R = RWS_Example_5_2()
    @test_logs (:warn,) KnuthBendix.knuthbendix(R, maxrules=200, implementation=:naive_kbs2)

    R = RWS_Example_5_3()
    test_kbs2_methods(R, (:naive_kbs2, :deletion), 6)
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs2)

    a,b = Word.([i] for i in 1:2)
    ε = one(a)
    @test Set(KnuthBendix.rules(rws)) == Set(KnuthBendix.Rule.([a^2=>ε, b^3=>ε,
        (b*a)^2=>a*b^2, a*b^2*a=>b*a*b, b^2*a*b^2=>a*b*a, (a*b)^2=>b^2*a]))

    R = RWS_Example_5_4()
    test_kbs2_methods(R, (:naive_kbs2, :deletion), 11)

    R = RWS_Example_5_5()
    test_kbs2_methods(R, (:naive_kbs2, :deletion), 18)

    R = RWS_Example_6_4()
    test_kbs2_methods(R, (:naive_kbs2, :deletion), 40)

    rws  = KnuthBendix.knuthbendix(R, maxrules=100, implementation=:naive_kbs2)
    rwsd = KnuthBendix.knuthbendix(R, maxrules=100, implementation=:deletion)

    @test Set(KnuthBendix.rules(rws)) == Set(KnuthBendix.rules(rwsd))

    w = Word([3, 3, 2, 2, 3, 3, 3, 1, 1, 1, 3, 1, 2, 3, 2, 3, 2, 3, 3, 3])

    @test KnuthBendix.rewrite_from_left(w, rws) == Word([1,3,1,2])

    R = RWS_Example_6_5()
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs2)
    @test _length(rws) == 40

    R = RWS_Closed_Orientable_Surface(3)
    rws = KnuthBendix.knuthbendix(R, implementation=:naive_kbs2)
    @test _length(rws) == 16
end

@testset "KBS-automata" begin
    for R in [
            RWS_Example_5_1(),
            # RWS_Example_5_2(), # non-confluent ℤ²
            RWS_Example_5_3(),
            RWS_Example_5_4(),
            RWS_Example_5_5(),
            # RWS_Example_6_4(),
            # RWS_Example_6_5(),
            RWS_Closed_Orientable_Surface(3),
        ]
        rws = KnuthBendix.knuthbendix2(R)
        R = KnuthBendix.knuthbendix2automaton!(R)
        @test _length(R) == _length(rws)
        @test Set(KnuthBendix.rules(rws)) == Set(KnuthBendix.rules(R))
    end
end


text:
type AccountingEngine{E <: AbstractMarketModelEvolver, P <: MarketModelMultiProduct}
  evolver::E
  product::P
  initialNumeraireValue::Float64
  numberProducts::Int
  numerairesHeld::Vector{Float64}
  numberCashFlowsThisStep::Vector{Int}
  cashFlowsGenerated::Vector{Vector{MarketModelCashFlow}}
  discounters::Vector{MarketModelDiscounter}
end

function AccountingEngine(evolver::AbstractMarketModelEvolver, product::MarketModelMultiProduct, initialNumeraireValue::Float64)
  product = clone(product)
  numProd = number_of_products(product)

  cashFlowsGenerated = Vector{MarketModelCashFlow}[[MarketModelCashFlow() for i = 1:max_number_of_cashflows_per_product_per_step(product)] for j = 1:numProd]

  cashFlowTimes = possible_cash_flow_times(product)
  rateTimes = get_evolution(product).rateTimes

  discounters = MarketModelDiscounter[MarketModelDiscounter(cashFlowTimes[i], rateTimes) for i in eachindex(cashFlowTimes)]

  return AccountingEngine(evolver, product, initialNumeraireValue, numProd, Vector{Float64}(numProd), Vector{Int}(numProd), cashFlowsGenerated, discounters)
end

function single_path_values!(ae::AccountingEngine, vals::Vector{Float64})
  fill!(ae.numerairesHeld, 0.0)
  weight = start_new_path!(ae.evolver)
  reset!(ae.product)

  principalInNumerairePortfolio = 1.0

  isDone = false

  while true
    thisStep = ae.evolver.currentStep
    weight *= advance_step!(ae.evolver)
    # println("forwards: ", ae.evolver.forwards)
    # println("logs: ", ae.evolver.logForwards)
    isDone = next_time_step!(ae.product, current_state(ae.evolver), ae.numberCashFlowsThisStep, ae.cashFlowsGenerated)

    numeraire = ae.evolver.numeraires[thisStep]

    # println("cfs: ", ae.cashFlowsGenerated)

    # for each product
    for i = 1:ae.numberProducts
      # and each cash flow
      cfs = ae.cashFlowsGenerated[i]

      for j = 1:ae.numberCashFlowsThisStep[i]
        # convert the cash flows to numeraires
        # this is done by calculating the number of numeraire bonds corresponding to such cash flow
        discounter = ae.discounters[cfs[j].timeIndex]
        bonds = cfs[j].amount * numeraire_bonds(discounter, current_state(ae.evolver), numeraire)

        # and adding the newly bought bonds to the number of numeraires held
        ae.numerairesHeld[i] += bonds / principalInNumerairePortfolio
      end
    end

    # error("die")

    if ~isDone
      # The numeraire might change between steps. This implies
      # that we might have to convert the numeraire bonds for
      # this step into a corresponding amount of numeraire
      # bonds for the next step. This can be done by changing
      # the principal of the numeraire and updating the number
      # of bonds in the numeraire portfolio accordingly.
      nextNumeraire = ae.evolver.numeraires[thisStep + 1]
      principalInNumerairePortfolio *= discount_ratio(current_state(ae.evolver), numeraire, nextNumeraire)
    end

    if isDone
      break
    end
  end

  for i in eachindex(ae.numerairesHeld)
    vals[i] = ae.numerairesHeld[i] * ae.initialNumeraireValue
  end

  return weight
end

function multiple_path_values!(ae::AccountingEngine, stats::GenericSequenceStats, numberOfPaths::Int)
  vals = Vector{Float64}(number_of_products(ae.product))
  QuantLib.Math.reset!(stats, length(vals), numberOfPaths)

  for i = 1:numberOfPaths
    weight = single_path_values!(ae, vals)
    # println("vals ", vals)
    # if i > 2
    #   error("DIE")
    # end
    add_sample!(stats, vals, i, weight)
  end

  return stats
end


text:
using Test
using ClimateMachine.BalanceLaws

import ClimateMachine.BalanceLaws: eq_tends, prognostic_vars

struct TestBL <: BalanceLaw end
struct X <: PrognosticVariable end
struct Y <: PrognosticVariable end
struct F1{PV} <: TendencyDef{Flux{FirstOrder}, PV} end
struct F2{PV} <: TendencyDef{Flux{SecondOrder}, PV} end
struct S{PV} <: TendencyDef{Source, PV} end

prognostic_vars(::TestBL) = (X(), Y())
eq_tends(::X, ::TestBL, ::Flux{FirstOrder}) = (F1{X}(),)
eq_tends(::Y, ::TestBL, ::Flux{FirstOrder}) = (F1{Y}(),)
eq_tends(::X, ::TestBL, ::Flux{SecondOrder}) = (F2{X}(),)
eq_tends(::Y, ::TestBL, ::Flux{SecondOrder}) = (F2{Y}(),)
eq_tends(::X, ::TestBL, ::Source) = (S{X}(),)
eq_tends(::Y, ::TestBL, ::Source) = (S{Y}(),)

@testset "BalanceLaws" begin
    bl = TestBL()
    @test prognostic_vars(bl) == (X(), Y())
    @test fluxes(bl, FirstOrder()) == (F1{X}(), F1{Y}())
    @test fluxes(bl, SecondOrder()) == (F2{X}(), F2{Y}())
    @test sources(bl) == (S{X}(), S{Y}())
    show_tendencies(bl)
    show_tendencies(bl; include_params = true)
end


text:
using Documenter, StackOverflow

makedocs(sitename="StackOverflow.jl")

deploydocs(
    repo = "github.com/logankilpatrick/StackOverflow.jl.git",
    target = "build"
)


text:
using Test, FinancialMonteCarlo,Statistics;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=[0.00,0.02];
T=1.0;
d=0.01;
r=0.02;
D=90.0;

Nsim=100000;
Nstep=30;
sigma=0.2;
mc=MonteCarloConfiguration(Nsim,Nstep);
toll=0.8

rfCurve=ZeroRate(r);

FwdData=Forward(T)
EUData=EuropeanOption(T,K)
Model_base=BrownianMotion(sigma,r-d-sigma*sigma*0.5,Underlying(0.0,0.0))
S=S0.*exp.(simulate(Model_base,mc,T));

Payoff=exp(-FinancialMonteCarlo.integral(r,T))*S[:,end];
@show Price=mean(Payoff)
Model=BlackScholesProcess(sigma,Underlying(S0,d));
@show FwdPrice=pricer(Model,rfCurve,mc,FwdData);


text:
using StregthMaterialsEdu
using Test

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


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

forwarddiff_dualized(::Type{TagType}, x::Real) where TagType = ForwardDiff.Dual{TagType}(x, one(x))

function forwarddiff_dualized(::Type{TagType}, x::NTuple{N,T}) where {TagType,N,T<:Real}
    ntuple(j -> ForwardDiff.Dual{TagType}(x[j], ntuple(i -> i == j ? one(x[j]) : zero(x[j]), Val(N))), Val(N))
end

# Equivalent to ForwardDiff internal function dualize(T, x):
forwarddiff_dualized(::Type{TagType}, x::SVector{N,T}) where {TagType,N,T<:Real} = SVector(forwarddiff_dualized(TagType, (x...,)))

# Equivalent to ForwardDiff internal function static_dual_eval(TagType, f, x) (for SVector):
function forwarddiff_eval(f::Base.Callable, x::Union{T, NTuple{N,T}, SVector{N,T}}) where {N,T<:Real}
    TagType = typeof(ForwardDiff.Tag(f, T))
    x_dual = forwarddiff_dualized(T, x)
    y_dual = f(x_dual)
end

forwarddiff_vjp(ΔΩ::Real, y_dual::ForwardDiff.Dual{TagType,T,1}) where {TagType,T<:Real} = ΔΩ * first(ForwardDiff.partials(y_dual))

function forwarddiff_vjp(ΔΩ::Union{NTuple{N,T},SVector{N,T}}, y_dual::NTuple{N,<:ForwardDiff.Dual}) where {N,T<:Real}
    (sum(map((ΔΩ_i, y_dual_i) -> ForwardDiff.partials(y_dual_i) * ΔΩ_i, ΔΩ, y_dual))...,)
end

# BAT.forwarddiff_vjp(ΔΩ, BAT.forwarddiff_eval(f, x)) == ForwardDiff.jacobian(f, x)' * ΔΩ (for SVector):
forwarddiff_vjp(ΔΩ::Union{NTuple{N,T},SVector{N,T}}, y_dual::SVector{N,<:ForwardDiff.Dual}) where {N,T<:Real} = SVector(forwarddiff_vjp((ΔΩ...,), (y_dual...,)))

# Return type of fwddiff_back (`SVector` or `NTuple`) currently depeds on type of `f(x)`, not type of `x`:
function forwarddiff_pullback(f::Base.Callable, x::Union{NTuple{N,T}, SVector{N,T}}) where {N,T<:Real}
    # Seems faster this way, according to benchmarking (benchmark artifact?):
    TagType = typeof(ForwardDiff.Tag(f, T))
    x_dual = forwarddiff_dualized(T, x)
    y_dual = f(x_dual)

    # Seems slower this way, for some reason:
    # y_dual = forwarddiff_eval(f, x)

    y = map(ForwardDiff.value, y_dual)
    fwddiff_back(ΔΩ) = forwarddiff_vjp(ΔΩ, y_dual)
    y, fwddiff_back
end

function forwarddiff_broadcast_pullback(fs, X::AbstractArray{<:Union{NTuple{N,T},SVector{N,T}}}) where {N,T<:Real}
    Y_dual = broadcast(forwarddiff_eval, fs, X)
    Y = broadcast(y_dual -> map(ForwardDiff.value, y_dual), Y_dual)
    fwddiff_back(ΔΩs) = broadcast(forwarddiff_vjp, ΔΩs, Y_dual)
    Y, fwddiff_back
end

#=
# May be faster for cheap `fs`, according to benchmarking, in spite of double evalutation (why?). Limited to SVector elements:
function forwarddiff_broadcast_pullback(fs, X::AbstractArray{<:SVector{N,T}}) where {N,T<:Real}
    Y = broadcast(fs, X)
    fwddiff_back(ΔΩs) = broadcast((f, x, ΔΩ) -> ForwardDiff.jacobian(f, x)' * ΔΩ, fs, X, ΔΩs)
    Y, fwddiff_back
end
=#


text:
using Compat

# First do the following from src/
#  include("tree.jl")
#  include("pvtree.jl")

# Make a tree like this:
# 0J
# 1C   2J   3C        4J    # children of 0
# 7C        5J   6J         # children of 1 and 3
# J = julia, C = C-code
# Node 7 is a garbage-collection
# With C == false, it should resolve to the following tree:
# 0J
# 2J   5J   6J   4J

node(ip::Integer) = PVTree.PVData(@compat(UInt(ip)), 1:0)

if VERSION < v"0.4.0-dev+980"
    isjl = Dict(uint(0:7),[true,false,true,false,true,true,true,false])
    isgc = Dict(uint(0:7),[falses(7);true])
else
    isjl = Dict(zip(@compat(map(UInt, 0:7)),[true,false,true,false,true,true,true,false]))
    isgc = Dict(zip(@compat(map(UInt, 0:7)),[falses(7);true]))
end

function buildraw()
    root = Tree.Node(node(0))
    c1 = Tree.addchild(root, node(1))
    Tree.addchild(root, node(2))
    c3 = Tree.addchild(root, node(3))
    Tree.addchild(root, node(4))
    Tree.addchild(c3, node(5))
    Tree.addchild(c3, node(6))
    Tree.addchild(c1, node(7))
    root
end

# Run it with C == false
root = buildraw()
PVTree.setstatus!(root, isgc)
PVTree.prunegraph!(root, isjl, (), (), (), [])

@assert root.data.status == 1
c = root.child
@assert c.data.ip == 2
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 5
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 6
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 4
@assert c.data.status == 0
@assert c == c.sibling

# Now do it again, with C == true
root = buildraw()
PVTree.setstatus!(root, isgc)
# PVTree.prunegraph!(root, isjl)

@assert root.data.status == 0
c = root.child
c1 = c
@assert c.data.ip == 1
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 2
@assert c.data.status == 0
c = c.sibling
c3 = c
@assert c.data.ip == 3
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 4
@assert c.data.status == 0
@assert c == c.sibling
c = c3.child
@assert c.data.ip == 5
@assert c.data.status == 0
c = c.sibling
@assert c.data.ip == 6
@assert c.data.status == 0
@assert c == c.sibling
c = c1.child
@assert c.data.ip == 7
@assert c.data.status == 1
@assert Tree.isleaf(c)
@assert c == c.sibling


text:
using RecurrenceAnalysis
using Base.Test

# Functions for dynamic systems
# Increment of x[ix] through RK4-integration
# dxdh is the array of derivatives, dt is the time increment
function dx_rk4(dxdt, x, ix, dt)
    x = collect(x)
    incr = zeros(length(x))
    k1 = dxdt[ix](x...)
    incr[ix] = dt*k1/2
    k2 = dxdt[ix]((x+incr)...)
    incr[ix] = dt*k2/2
    k3 = dxdt[ix]((x+incr)...)
    incr[ix] = dt*k3
    k4 = dxdt[ix]((x+incr)...)
    dt/6*(k1 + 2k2 + 2k3 + k4)
end
# Integrate dynamical system through n intervals of length dt
function dynamical_system(x0, dxdt, dt, n)
    m = length(x0)
    x = zeros(n,m)
    x[1,:] = collect(x0)
    for t=1:n-1
        dx = zeros(m)
        for ix=1:m
            dx[ix] = dx_rk4(dxdt, x[t,:], ix, dt)
        end
        x[t+1,:] = x[t,:][:] + dx
    end
    x
end

# Test with Lorenz system
lorenz_eq(sigma, rho, b) = (
  (x,y,z) -> sigma*(y - x),
  (x,y,z) -> x*(rho - z) - y,
  (x,y,z) -> x*y - b*z)

sigma=10.; rho=28.; b=8/3;
x0 = ones(3)
lorenz_data = dynamical_system(x0, lorenz_eq(sigma,rho,b), .01, 1000)
x = lorenz_data[501:2:end,1]
# Look for optimal delay
ami_def = ami(x, (1,12))
@test findmin(ami_def)[2] == 9
ami_fd = ami(x, (1,12), "FD")
@test findmin(ami_fd)[2] == 10
ami_15 = ami(x, (1,12), 15)
gmi_10 = gmi(x, (1,12), 0.1)
# Look for optimal embedding dimension
fnnval  =  fnn(x, (1,5), 8, (15, 2))
e1,e2   = afnn(x, (1,5), 8)
ffnnval = ffnn(x, (1,5), 8)
# Look for optimal threshold
dd, rr = sorteddistances(x, theiler=1)
# Distance and recurrence matrices
xe = embed(x, 3, 8)
dmat = distancematrix(xe)
dmat_euc = distancematrix(xe,"euclidean")
@test all(dmat .<= dmat_euc)
rmat = recurrencematrix(xe, 0.05)
crmat = crossrecurrencematrix(x, lorenz_data[:,3], 0.05)
jrmat = jointrecurrencematrix(x, lorenz_data[:,3], 0.05)
# RQA
rr = recurrencerate(rmat)
@test recurrencerate(rmat,theiler=1) < rr
detm = determinism(rmat)
@test determinism(rmat, theiler=5, lmin=5) < detm
lmean = avgdiag(rmat)
@test avgdiag(rmat, lmin=5) > lmean
@test maxdiag(rmat) > lmean
ent = entropy(rmat)
tnd = trend(rmat, border=20, theiler=5)
lam = laminarity(rmat)
tt = trappingtime(rmat, lmin=5)
@test maxvert(rmat) >= tt


text:
"""
Provides a definition of a FULL method for constructing Koza trees.
"""
type KozaFullBuilderDefinition <: KozaBuilderDefinition
  """
  The minimum depth of trees created according to this method.
  """
  min_depth::Int

  """
  The maximum depth of trees created according to this method.
  """
  max_depth::Int

  KozaFullBuilderDefinition() = new(1, 8)
  KozaFullBuilderDefinition(min::Int, max::Int) = new(min, max)
end

"""
Implements the Koza FULL initialisation method.
"""
type KozaFullBuilder{T <: KozaTree} <: KozaBuilder
  """
  The minimum depth of trees created according to this method.
  """
  min_depth::Int

  """
  The maximum depth of trees created according to this method.
  """
  max_depth::Int

  KozaFullBuilder() = new(1, 8)
  KozaFullBuilder(min::Int, max::Int) = new(min, max)
end

"""
Composes a Koza Full builder, for a given sub-type of Koza Tree, using a
provided definition.
"""
compose!(d::KozaFullBuilderDefinition, t::Type) =
  KozaFullBuilder{t}(d.min_depth, d.max_depth)

"""
Implements the Koza FULL initialisation method, where some `d` is selected
between a specified minimum and maximum depth, and used to generate a full tree
with a depth equal to `d`.

**Properties:**

* `min::Int`, the minimum tree depth.
* `max::Int`, the maximum tree depth.
"""
full() = KozaFullBuilderDefinition()
full(min::Int, max::Int) = KozaFullBuilderDefinition(min, max)
function full(f::Function)
  def = KozaFullBuilderDefinition()
  f(def)
  def
end

"""
Constructs a Koza Tree belonging to a given representation using the FULL
method, according to a set of provided specifics.
"""
function build{T}(b::KozaFullBuilder{T}, r::KozaTreeRepresentation)
  t = T()
  t.root = build(KozaFullBuilder, r, t, 1, rand(b.min_depth:b.max_depth))
  return t
end

build(::Type{KozaFullBuilder},
  r::KozaTreeRepresentation,
  p::KozaParent,
  d::Integer,
  md::Integer
) =
  if d == md
    return sample_terminal(r, p)
  else
    n = sample_non_terminal(r, p)
    for i in 1:arity(n)
      n.children[i] = build(KozaFullBuilder, r, n, d+1, md)
    end
    return n
  end


text:
@testset "map" begin
    @inline foo(x, y) = exp(x) - sin(y)
    for T ∈ (Float32,Float64)
        @show T, @__LINE__
        for N ∈ [ 3, 371 ]
            a = rand(T, N); b = rand(T, N);
            c1 = map(foo, a, b);
            c2 = vmap(foo, a, b);
            @test c1 ≈ c2
            c2 = vmapnt(foo, a, b);
            @test c1 ≈ c2
            c2 = vmapntt(foo, a, b);
            @test c1 ≈ c2
            fill!(c2, NaN); @views vmapnt!(foo, c2[2:end], a[2:end], b[2:end]);
            @test @views c1[2:end] ≈ c2[2:end]
            fill!(c2, NaN); @views vmapntt!(foo, c2[2:end], a[2:end], b[2:end]);
            @test @views c1[2:end] ≈ c2[2:end]
        end

        c = rand(T,100); x = rand(T,10^4); y1 = similar(x); y2 = similar(x);
        map!(xᵢ -> clenshaw(xᵢ, c), y1, x)
        vmap!(xᵢ -> clenshaw(xᵢ, c), y2, x)
        @test y1 ≈ y2

    end
end


text:
"""
    Figure(workstation::Tuple{Float64, Float64}, plots::Vector{PlotObject})

Return a new figure, defined by:

* **`workstation`**: a `Tuple{Float64, Float64}` with the width and height of the
    overall plot container (workstation), in pixels.
* **`plots`**: a vector of [`PlotObject`](@ref) elements, which contain the information of
    the individual plots included in the figure.
"""
struct Figure
    workstation::Tuple{Float64, Float64}
    plots::Vector{PlotObject}
end

const EMPTYFIGURE = Figure((0.0, 0.0), [PlotObject()])
const CURRENTFIGURE = Ref(EMPTYFIGURE)

"""
    gcf()
    gcf(fig::Figure)

Get the global current figure, or set it to be `fig`.
"""
gcf() = (CURRENTFIGURE[] == EMPTYFIGURE) ? Figure() : CURRENTFIGURE[]
gcf(fig::Figure) = (CURRENTFIGURE[] = fig)

"""
    Figure([figsize::Tuple{Float64, Float64}, units::String])

Create a new figure of a given size.

The figure size is defined by `figsize` (a 2-tuple with the target width and height),
and `units` (a string with the abbreviation of the units: "px" for pixels,
"in" for inches, "cm" for centimeters or "m" for meters).
The default dimensions used by `Figure()` with no arguments is 600×450 pixels
— or a proportional increased size, if the detected display resolution is high.
This constructor also sets the "current figure" to the one that has just been
created..
"""
function Figure(figsize=(600,450), units::String="px")
    # Calculate the display resolution (dpi)
    mwidth, mheight, width, height = GR.inqdspsize()
    dpi = width / mwidth * 0.0254
    # Transform figsize to pixels
    sizepx = Tuple(float(s) for s in figsize)
    if units == "px" # keep it as is
    elseif units == "in"
        sizepx = figsize .* dpi
    elseif units == "cm"
        sizepx = figsize .* (dpi / 2.54)
    elseif units == "m"
        sizepx = figsize .* (dpi / 0.0254)
    end
    # increase dimensions if dpi is too big (the figure would be tiny)
    w, h = (dpi > 200) ? Tuple(s * dpi / 100 for s in sizepx) : sizepx
    # Workstation size
    wssize = mwidth / width * w
    ratio = (w > h) ? (1.0, float(h)/w) : (float(w)/h, 1.0)
    workstation = (wssize * ratio[1], wssize * ratio[2])
    plots = [PlotObject()]
    CURRENTFIGURE[] = Figure(workstation, plots)
end

function Base.Multimedia.display(fig::Figure)
    output = draw(fig)
    output isa Nothing && return nothing
    display(output)
    return nothing
end

# Caution! This depends on GR internals
# Base.showable(::MIME"image/svg+xml", ::Figure) = GR.mime_type[] == "svg"
# Base.showable(::MIME"image/png", ::Figure) = GR.mime_type[] == "png"
Base.showable(::MIME"text/html", ::Figure) = GR.mime_type[] ∈ ("mov", "mp4", "webm")

# Base.show(io::IO, mime::M, fig::Figure) where {
#     M <: Union{MIME"image/svg+xml", MIME"image/png", MIME"text/html"}
#     } = show(io, mime, draw(fig))

Base.show(io::IO, mime::MIME"text/html", fig::Figure) = show(io, mime, draw(fig))

for (mime, fmt) in (
    "image/png" => "png",
    "image/svg+xml" => "svg"
)

    @eval function Base.show(io::IO, mime::MIME{Symbol($mime)}, fig::Figure)
        gr_mime = GR.isinline() ? string(GR.mime_type[]) : ""
        GR.inline($fmt)
        try
            show(io, mime, draw(fig))
        catch err
            throw(err)
        finally
            resetmime(gr_mime)
        end
    end
end

"""
    currentplot([fig::Figure])

Get the "current plot", i.e. the target of the next plotting operations in the
(optionally) given figure `fig`.

If no figure is given, the "current figure" is used (cf. [`gcf`](@ref)).
"""
currentplot(f::Figure=gcf()) = f.plots[end]

# Normalized width and height of a figure'w workstation
wswindow(f::Figure) = f.workstation ./ maximum(f.workstation)

## Subplots

function subplot!(f::Figure, nr, nc, p, replace=true)
    xmin, xmax, ymin, ymax = 1.0, 0.0, 1.0, 0.0
    for i in collect(p)
        r = nr - div(i-1, nc)
        c = (i-1) % nc + 1
        xmin = min(xmin, (c-1)/nc)
        xmax = max(xmax, c/nc)
        ymin = min(ymin, (r-1)/nr)
        ymax = max(ymax, r/nr)
    end
    coords = [xmin, xmax, ymin, ymax]
    po = PlotObject(; subplot = coords)
    if replace
        po = replaceplot!(f.plots, po)
    end
    push!(f.plots, po)
    return po
end

"""
    subplot(num_rows, num_columns, indices[, replace])

Set a subplot in the current figure.

By default, the current plot covers the whole window. To display more
than one plot, the window can be split into a number of rows and columns,
with each plot covering one or more cells in the resulting grid.

Subplot indices are one-based and start at the upper left corner, with a
new row starting after every `num_cols` subplots.

The arguments `num_rows` and `num_cols` indicate the number of rows and columns
of the grid of plots into which the figure is meant to be divided, and `indices`
is an integer or a collection of integers that identify a group of cells in that
grid. This function returns a plot with the minimum size that spans over all
those cells, and appends it to the array of plots of the figure,
such that it becomes its current plot.

If the viewport of the new subplot coincides with the viewport of an existing
plot, by default the older one is moved to the first plane, and taken as the
"current plot"; but if there is only a partial overlap between the new subplot
and other plots, the overlapping plots are removed.
To override this behavior and keep all previous plots without changes,
set the optional argument `replace` to false.

# Examples

```julia
$(_example("subplot"))
```
"""
subplot(args...) = subplot!(gcf(), args...)

function replaceplot!(plotcollection, p)
    # If subplot is not specified in p, empty the whole collection
    if !haskey(p.attributes, :subplot)
        empty!(plotcollection)
    else
        coords = p.attributes[:subplot]
        i = 1
        while i ≤ length(plotcollection)
            # Check intersection
            coords_i = get(plotcollection[i].attributes, :subplot, [0.0, 1.0, 0.0, 1.0])
            bottomleft = (max(coords[1], coords_i[1]), max(coords[3], coords_i[3]))
            topright = (min(coords[2], coords_i[2]), min(coords[4], coords_i[4]))
            if all(bottomleft .< topright)
                (coords ≈ coords_i) && (p = plotcollection[i])
                deleteat!(plotcollection, i)
            else
                i += 1
            end
        end
    end
    return p
end

function draw(f::Figure)
    GR.clearws()
    w, h = f.workstation
    GR.setwsviewport(0.0, w, 0.0, h)
    ratio_w, ratio_h = wswindow(f)
    GR.setwswindow(0.0, ratio_w, 0.0, ratio_h)
    drawn = false
    for p in f.plots
        drawn = draw(p) || drawn
    end
    GR.updatews()
    if GR.isinline() && drawn
        return GR.show()
    else
        return
    end
end


text:
function make_count_matrix(counts)
    idx = findall(counts .> 0)
    countmat = Array{Int}(undef, length(idx), 2)
    countmat[:, 1] = idx
    countmat[:, 2] = counts[idx]
    return countmat
end

function format_counts_lda(countsdf::DataFrame, cols::Vector{Symbol})
    counts = Matrix{Int}[]
    for col in cols
        colcounts = convert(Array, countsdf[!, col])
        countmat = make_count_matrix(colcounts)

        push!(counts, countmat)
    end
    return counts
end

function format_counts_ctm(countsdf::DataFrame, cols::Vector{Symbol})
    return format_counts_mmctm([countsdf], cols)
end

function format_counts_mmctm(countdfs::Vector{DataFrame}, cols::Vector{Symbol})
    counts = Vector{Matrix{Int}}[]
    for col in cols
        doc_counts = Matrix{Int}[]
        for df in countdfs
            push!(doc_counts, make_count_matrix(convert(Array, df[!, col])))
        end

        push!(counts, doc_counts)
    end

    return counts
end


text:
#####
##### Upwind-biased 1rd-order advection scheme
#####

"""
    struct UpwindBiasedFirstOrder <: AbstractUpwindBiasedAdvectionScheme{1}

Upwind-biased first-order advection scheme.
"""
struct UpwindBiasedFirstOrder <: AbstractUpwindBiasedAdvectionScheme{1} end

const U1 = UpwindBiasedFirstOrder

@inline boundary_buffer(::U1) = 1

@inline symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, ::U1, c) = ℑxᶠᵃᵃ(i, j, k, grid, c)
@inline symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, ::U1, c) = ℑyᵃᶠᵃ(i, j, k, grid, c)
@inline symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, ::U1, c) = ℑzᵃᵃᶠ(i, j, k, grid, c)

@inline symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, ::U1, u) = ℑxᶜᵃᵃ(i, j, k, grid, u)
@inline symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, ::U1, v) = ℑyᵃᶜᵃ(i, j, k, grid, v)
@inline symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, ::U1, w) = ℑzᵃᵃᶜ(i, j, k, grid, w)

@inline left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, ::U1, c) = @inbounds c[i-1, j, k]
@inline left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, ::U1, c) = @inbounds c[i, j-1, k]
@inline left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, ::U1, c) = @inbounds c[i, j, k-1]

@inline left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::U1, u) = left_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, u)
@inline left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::U1, v) = left_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, v)
@inline left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::U1, w) = left_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, w)

@inline right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, ::U1, c) = @inbounds c[i, j, k]
@inline right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, ::U1, c) = @inbounds c[i, j, k]
@inline right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, ::U1, c) = @inbounds c[i, j, k]

@inline right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::U1, u) = right_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, u)
@inline right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::U1, v) = right_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, v)
@inline right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::U1, w) = right_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, w)


text:
@enum TokenIDS begin
   IDENTIFIER     = 0
   ASSIGN         = 1
   INT_NUMBER     = 2
   FLOAT_NUMBER   = 3
   OPERATOR       = 4
   RESERVED_WORD  = 5
   PUNCTUATION    = 6
   CHAR           = 7
   STRING         = 8
   TYPE           = 9
   INVALID        = 10
   WHITESPACE     = 11
   CFLUX          = 12
   BOOL           = 13
   COMMENT        = 14
   EOF            = -1
end # enum

mutable struct Token
    id   :: TokenIDS
    text :: String
    span :: Tuple{Integer, Integer}
    line :: Integer
    col  :: Integer
    function Token(id=INVALID, txt="", span=(0, 0), line=1, col=1)
        new(id, txt, span, line, col)
    end # function
end # struct

mutable struct Tokens
   tokens :: Array{Token}
   pos    :: Integer
   total  :: Integer
   function Tokens(tks::Array{Token})
      new(tks, 0, length(tks))
   end
end

function next!(t::Tokens)
   if t.pos == t.total
      error("EOF reached and parser not finished", previous(t))
   end

   t.pos += 1
   v = t.tokens[ t.pos ]

   if v.id == INVALID
      error("Invalid token $(v.text)", v)
   end
   return v
end

function reset!(t::Tokens)
   t.pos = 0
end

next(t::Tokens)      = t.tokens[(t.pos + 1)]
previous(t::Tokens)  = t.tokens[(t.pos - 1)]
current(t::Tokens)   = t.tokens[t.pos]
roll_back(t::Tokens) = t.pos -= 1

# Base.show is invoked in Atom and @show
Base.show(io::IO, t::Token) = println(io, "`$(t.text)` is $(t.id) start: $(t.span[1]) end: $(t.span[2]) at $(t.line):$(t.col)")
# Base.println is invoked in println call
Base.println(t::Token) = println("`$(t.text)` is $(t.id) start: $(t.span[1]) end: $(t.span[2]) at $(t.line):$(t.col)")

max_stars_repo_name: GabrielTrettel/Portran

max_issues_repo_name: GabrielTrettel/Portran

max_forks_repo_name: GabrielTrettel/Portran


text:
using Revise

using BNPRegressionGGA2021
using Distributions
using Gadfly
using Random
using Colors

const BNP = BNPRegressionGGA2021

#=================================#
# Example 1 - 1 trivial predictor #
#=================================#

function simulate_sample(N0, N1)
    dy = MixtureModel(Normal, [(1.0, 1.0), (-1.0, 0.25)], [0.5, 0.5])
    X0 = ones(N0, 1)
    X1 = ones(N1, 1)
    y0 = rand(dy, N0)
    y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect
    return dy, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy, y0, X0, y1, X1 = simulate_sample(N0, N1);
m = BNP.DGSBPNormal(; y0, X0, y1, X1);
chainf, chainβ, chaing = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);

# With Gadfly
f1 = mean(chainf);
f0 = pdf.(dy, y1);
p = plot(
    layer(x = y1, y = f1, Geom.line, color = ["bnp"]),
    layer(x = y0, Geom.density, color = ["kden"]),
    layer(x = y0, Geom.histogram(density = true), color = ["hist"]),
    style(major_label_font="Arial", minor_label_font="Arial"),
)

# #========================================================#
# # Example 2 - 1 discrete predictor (simple distribution) #
# #========================================================#

# # Simulate a sample
# function simulate_sample(N0, N1)
#     x0 = [-0.5, 0.5]
#     X0 = [ones(N0) rand(x0, N0)]
#     y0 = X0 * ones(2) + randn(N0)
#     y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect |> x -> repeat(x, 2)
#     x1 = [x0[1] * ones(N1); x0[2] * ones(N1)]
#     X1 = [ones(2 * N1) x1]
#     return y0, X0, y1, X1
# end

# Random.seed!(1);
# N0, N1 = 500, 50
# y0, X0, y1, X1 = simulate_sample(N0, N1);
# m = BNP.DGSBPNormal(; y0, X0, y1, X1);
# chainf, chainβ = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);

# fb = mean(chainf);
# plot(x = y0, color = X0[:, 2], Geom.density)
# plot(x = y1, y = fb, color = X1[:, 2], Geom.line)
# mean(chainβ)

# #=================================================#
# # Example 3 - 2 predictors (mixture distribution) #
# #=================================================#

# # Simulate a sample
# function simulate_sample(N0, N1)
#     dϵ = MixtureModel(Normal, [(1.0, 1.0), (-1.0, 0.25)], [0.5, 0.5])
#     x0 = [-0.5, 0.5]
#     X0 = [ones(N0) rand(x0, N0, 2)]
#     y0 = X0 * [ones(2); 0] + rand(dϵ, N0)
#     y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect |> x -> repeat(x, 4)
#     x1 = kron([x0[1], x0[2], x0[1], x0[2]], ones(N1))
#     x2 = kron([x0[1], x0[1], x0[2], x0[2]], ones(N1))
#     X1 = [ones(4 * N1) x1 x2]
#     return y0, X0, y1, X1
# end

# Random.seed!(1);
# N0, N1 = 500, 50
# y0, X0, y1, X1 = simulate_sample(N0, N1);
# m = BNP.DGSBPNormal(; y0, X0, y1, X1);
# chainf, chainβ = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);

# fb = mean(chainf);
# X0concat = @. string(X0[:, 2]) * " " * string(X0[:, 3]);
# X1concat = @. string(X1[:, 2]) * " " * string(X1[:, 3]);
# plot(x = y0, color = X0concat, Geom.density)
# plot(x = y1, y = fb, color = X1concat, Geom.line)
# mean(chainβ)

# #=================================================#
# # Example 4 - 3 predictors (mixture distribution) #
# #=================================================#

# # Simulate a sample
# function simulate_sample(N0, N1)
#     dϵ = MixtureModel(Normal, [(1.0, 1.0), (-1.0, 0.25)], [0.5, 0.5])
#     x0 = [-0.5, 0.5]
#     X0 = [ones(N0) rand(x0, N0, 3)]
#     y0 = X0 * [1, 0, 1, 0] + rand(dϵ, N0)
#     y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect |> x -> repeat(x, 8)
#     x1 = kron([x0[1], x0[1], x0[1], x0[1], x0[2], x0[2], x0[2], x0[2]], ones(N1))
#     x2 = kron([x0[1], x0[1], x0[2], x0[2], x0[1], x0[1], x0[2], x0[2]], ones(N1))
#     x3 = kron([x0[1], x0[2], x0[1], x0[2], x0[1], x0[2], x0[1], x0[2]], ones(N1))
#     X1 = [ones(8 * N1) x1 x2 x3]
#     return y0, X0, y1, X1
# end

# Random.seed!(1);
# N0, N1 = 500, 50
# y0, X0, y1, X1 = simulate_sample(N0, N1);
# m = BNP.DGSBPNormal(; y0, X0, y1, X1);
# chainf, chainβ = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);

# fb = mean(chainf);
# X0concat = @. string(X0[:, 2]) * " " * string(X0[:, 3]) * " " * string(X0[:, 4]);
# X1concat = @. string(X1[:, 2]) * " " * string(X1[:, 3]) * " " * string(X1[:, 4]);
# plot(x = y0, color = X0concat, Geom.density)
# plot(x = y1, y = fb, color = X1concat, Geom.line)
# mean(chainβ)

# #==================================================#
# # Example 5 - 10 predictors (mixture distribution) #
# #==================================================#

# # Simulate a sample
# function simulate_sample(N0, N1)
#     dϵ = MixtureModel(Normal, [(1.0, 1.0), (-1.0, 0.25)], [0.5, 0.5])
#     x0 = [-0.5, 0.5]
#     X0 = [ones(N0) rand(x0, N0, 10)]
#     y0 = X0 * [1; zeros(9); 1] + rand(dϵ, N0)
#     y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect |> x -> repeat(x, 2)
#     x1 = kron([x0[1], x0[2]], ones(N1))
#     X1 = [ones(2 * N1) x1 ones(2 * N1, 9)]
#     return y0, X0, y1, X1
# end

# Random.seed!(1);
# N0, N1 = 500, 50
# y0, X0, y1, X1 = simulate_sample(N0, N1);
# m = BNP.DGSBPNormal(; y0, X0, y1, X1);
# chainf, chainβ = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);
# mean([chainβ[i] .== zeros(11) for i in 1:length(chainβ)])


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((false, 1, 2, true, true),),
    ((false, 1, 2, false, true),),
    ((false, 1, 2, false, false),),
    ((false, 3, 1, false, false),),
    ((true, :abs1d, 1, true, true),),
    ((true, :absball2d, 1, false, true),),
    ((false, :motzkin, 3, true, true),),
    ((false, :motzkin_ball, 3, false, false),),
    ((false, :motzkin_ellipsoid, 3, false, true),),
    ]
insts["fast"] = [
    ((false, 1, 30, true, true),),
    ((false, 1, 30, false, true),),
    ((false, 1, 30, false, false),),
    ((false, 2, 8, true, true),),
    ((false, 3, 6, true, true),),
    ((false, 5, 3, true, true),),
    ((false, 10, 1, true, true),),
    ((false, 10, 1, false, true),),
    ((false, 10, 1, false, false),),
    ((false, 4, 4, true, true),),
    ((false, 4, 4, false, true),),
    ((false, 4, 4, false, false),),
    ((true, :abs1d, 3, true, true),),
    ((true, :absunit1d, 1, true, true),),
    ((true, :absunit1d, 3, true, true),),
    ((true, :negabsunit1d, 2, true, true),),
    ((true, :absball2d, 1, true, true),),
    ((true, :absbox2d, 2, true, true),),
    ((true, :negabsbox2d, 1, true, true),),
    ((true, :denseunit1d, 2, true, true),),
    ((true, :negabsunit1d, 2, false, true),),
    ((true, :absball2d, 1, false, true),),
    ((true, :negabsbox2d, 1, false, true),),
    ((true, :denseunit1d, 2, false, true),),
    ((false, :butcher, 2, true, true),),
    ((false, :caprasse, 4, true, true),),
    ((false, :goldsteinprice, 7, true, true),),
    ((false, :goldsteinprice_ball, 6, true, true),),
    ((false, :goldsteinprice_ellipsoid, 7, true, true),),
    ((false, :heart, 2, true, true),),
    ((false, :lotkavolterra, 3, true, true),),
    ((false, :magnetism7, 2, true, true),),
    ((false, :magnetism7_ball, 2, true, true),),
    ((false, :motzkin, 3, true, true),),
    ((false, :motzkin_ball, 3, true, true),),
    ((false, :motzkin_ellipsoid, 3, true, true),),
    ((false, :reactiondiffusion, 4, true, true),),
    ((false, :robinson, 8, true, true),),
    ((false, :robinson_ball, 8, true, true),),
    ((false, :rosenbrock, 5, true, true),),
    ((false, :rosenbrock_ball, 5, true, true),),
    ((false, :schwefel, 2, true, true),),
    ((false, :schwefel_ball, 2, true, true),),
    ((false, :lotkavolterra, 3, false, true),),
    ((false, :motzkin, 3, false, true),),
    ((false, :motzkin_ball, 3, false, true),),
    ((false, :schwefel, 2, false, true),),
    ((false, :lotkavolterra, 3, false, false),),
    ((false, :motzkin, 3, false, false),),
    ((false, :motzkin_ball, 3, false, false),),
    ]
return (PolyMinNative, insts)


text:
using LightGraphs

function build_graph(file::AbstractString)
    objects = Set{String}()
    for line in eachline(file), object in split(line, ')')
        push!(objects, object)
    end
    vertex_map = Dict(object => i for (i, object) in enumerate(objects))
    graph = SimpleDiGraph(length(objects))
    for line in eachline(file)
        a, b = split(line, ')')
        add_edge!(graph, vertex_map[a], vertex_map[b])
    end
    (graph, vertex_map)
end

function main(input)
    directed_graph, vertex_map = build_graph(input)
    i = first(first(attracting_components(reverse(directed_graph))))
    println("Part 1 answer: ", sum(gdistances(directed_graph, i)))

    graph = SimpleGraph(directed_graph)
    you = vertex_map["YOU"]
    you_parent = first(outneighbors(graph, you))
    santa = vertex_map["SAN"]
    santa_parent = first(outneighbors(graph, santa))
    println("Part 2 answer: ", length(a_star(graph, you_parent, santa_parent)))
end

main("input")

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
# TODO: Implement the following as a performant (when the number of columns is
# large) alternative to a zipped iterator
# """
# """
# immutable TableRowIterator{T}
#     cols::T
# end
#
# Base.start(itr::TableRowIterator) = 1
# @generated function Base.next{T}(itr::TableRowIterator{T}, st)
#     ncols = length(T.parameters)
#     tuple_ex = Expr(:tuple)
#     for j in 1:ncols
#         push!(tuple_ex.args, :( itr.cols[$j][st] ))
#     end
#     return quote
#         return ($tuple_ex, st+1)
#     end
# end
# Base.done(itr::TableRowIterator, st) =
#     length(itr.cols) > 0 ? done(itr.cols[1], st) : true
#
# AbstractTables.eachrow(tbl::Table) = TableRowIterator(tuple(columns(tbl)...))
# function AbstractTables.eachrow(tbl::Table, flds...)
#     cols = [ tbl[fld] for fld in flds ]
#     TableRowIterator(tuple(cols...))
# end

"""
"Enumerate" the columns of an AbstractTable by field (column name).

Arguments:

* tbl::AbstractTable

Returns:

* cols::Base.Zip2{Array{Symbol,1},Array{Any,1}}: An iterator over
`(field, column)` pairs from `tbl`.
"""
eachcol(tbl::AbstractTable) = zip(fields(tbl), columns(tbl))

"""
    `eachrow(tbl[, flds...])`

Return an iterator over the rows (materialized as `Tuple`s) of `tbl`. If `flds`
are specified, the tuples returned by the row iterator will contain data from
only the respective columns.
"""
eachrow(tbl::AbstractTable) = zip(columns(tbl)...)
function eachrow(tbl::AbstractTable, fields...)
    cols = [ tbl[field] for field in fields ]
    return zip(cols...)
end


text:
include("MetaDict.jl")

export Image, ImageHeader, ImageFlags

@MRD.enm ImageType::UInt16 magnitude = 1 phase = 2 real = 3 imag = 4 complex = 5
@MRD.enm TypeIndex::UInt16 ushort = 1 short = 2 uint = 3 int = 4 float = 5 double = 6 cxfloat =
    7 cxdouble = 8

datatype_to_type = Base.Dict{TypeIndex.Enm,Type}([
    (TypeIndex.ushort, UInt16),
    (TypeIndex.short, Int16),
    (TypeIndex.uint, UInt32),
    (TypeIndex.int, Int32),
    (TypeIndex.float, Float32),
    (TypeIndex.double, Float64),
    (TypeIndex.cxfloat, ComplexF32),
    (TypeIndex.cxdouble, ComplexF64),
])
type_to_datatype = Base.Dict(value => key for (key, value) in datatype_to_type)
is_mrd_datatype(X) = X ∈ keys(type_to_datatype)

module ImageFlags
import FlagSets
FlagSets.@flagset Flags {_, UInt64} begin
    :is_navigation_data
    72057594037927936 --> :user1
    :user2
    :user3
    :user4
    :user5
    :user6
    :user7
    :user8
end
end

@with_kw struct ImageHeader
    version::UInt16 = 1
    flags::ImageFlags.Flags = ImageFlags.Flags(0)
    measurement_uid::UInt32 = 0
    field_of_view::NTuple{3,Float32} = (0, 0, 0)
    position::NTuple{3,Float32} = (0, 0, 0)
    read_dir::NTuple{3,Float32} = (1, 0, 0)
    phase_dir::NTuple{3,Float32} = (0, 1, 0)
    slice_dir::NTuple{3,Float32} = (0, 0, 1)
    patient_table_position::NTuple{3,Float32} = (0, 0, 0)
    average::UInt16 = 0
    slice::UInt16 = 0
    contrast::UInt16 = 0
    phase::UInt16 = 0
    repetition::UInt16 = 0
    set::UInt16 = 0
    acquisition_time_stamp::UInt32 = 0
    physiology_time_stamp::NTuple{3,UInt32} = (0, 0, 0)
    image_type::ImageType.Enm = ImageType.magnitude
    image_index::UInt16 = 0
    image_series_index::UInt16 = 0
    user_int::NTuple{8,Int32} = (0, 0, 0, 0, 0, 0, 0, 0)
    user_float::NTuple{8,Float32} = (0, 0, 0, 0, 0, 0, 0, 0)
end

struct Image
    header::ImageHeader
    data::Array{T,4} where {T<:Number}
    meta::MetaDict
    function Image(header::ImageHeader, data::Array, meta = MetaDict())
        is_mrd_datatype(eltype(data)) || error("Only MRD standard types supported")
        ndims(data) <= 4 || error("Images have a maximum of 4 dimensions")
        data = reshape(data, Val(4))
        new(header, data, meta)
    end

end

function Base.:(==)(img1::Image,img2::Image)
    return img1.header == img2.header && img1.data == img2.data && img1.meta == img2.meta
end

struct RawImageHeader

    version::UInt16
    data_type::TypeIndex.Enm
    flags::ImageFlags.Flags
    measurement_uid::UInt32
    matrix_size::NTuple{3,UInt16}
    field_of_view::NTuple{3,Float32}
    channels::UInt16
    position::NTuple{3,Float32}
    read_dir::NTuple{3,Float32}
    phase_dir::NTuple{3,Float32}
    slice_dir::NTuple{3,Float32}
    patient_table_position::NTuple{3,Float32}
    average::UInt16
    slice::UInt16
    contrast::UInt16
    phase::UInt16
    repetition::UInt16
    set::UInt16
    acquisition_time_stamp::UInt32
    physiology_time_stamp::NTuple{3,UInt32}
    image_type::ImageType.Enm
    image_index::UInt16
    image_series_index::UInt16
    user_int::NTuple{8,Int32}
    user_float::NTuple{8,Float32}
    attribute_string_len::UInt32
end

function RawImageHeader(image::Image)
    channels = UInt16(size(image.data)[4])
    matrix_size = Tuple{UInt16,UInt16,UInt16}(size(image.data)[1:3])
    data_type = type_to_datatype[eltype(image.data)]
    hdr = image.header
    attribute_string_len = ncodeunits(to_xml_string(image.meta))

    RawImageHeader(
        hdr.version,
        data_type,
        hdr.flags,
        hdr.measurement_uid,
        matrix_size,
        hdr.field_of_view,
        channels,
        hdr.position,
        hdr.read_dir,
        hdr.phase_dir,
        hdr.slice_dir,
        hdr.patient_table_position,
        hdr.average,
        hdr.slice,
        hdr.contrast,
        hdr.phase,
        hdr.repetition,
        hdr.set,
        hdr.acquisition_time_stamp,
        hdr.physiology_time_stamp,
        hdr.image_type,
        hdr.image_index,
        hdr.image_series_index,
        hdr.user_int,
        hdr.user_float,
        attribute_string_len
    )

end

@generated function ImageHeader(raw::RawImageHeader)
    names = fieldnames(ImageHeader)
    args = names .|> s -> Expr(:call,:getfield,:raw,QuoteNode(s))
    :(ImageHeader($(args...)))
end

function ImageHeader(acq::AcquisitionHeader; kwargs...)
    names = (:version, :measurement_uid, :position, :read_dir, :phase_dir, :slice_dir, :patient_table_position, :acquisition_time_stamp, :physiology_time_stamp, :user_int, :user_float )

    fields_acq = LittleDict(n => getfield(acq,n) for n in names)

    names_idx = (:average,:slice, :contrast,:phase,:repetition, :set)
    fields_idx = LittleDict(n => getfield(acq.idx,n) for n in names_idx)

    fields = merge(fields_acq,fields_idx,kwargs)

    return ImageHeader(;fields...)
end

write(io::IO, header::RawImageHeader) = fieldnames(RawImageHeader) .|> getfield $ header .|> MRD.unsafe_write $ io

function write(io::IO, img::Image)
    write(io, RawImageHeader(img))
    write(io,img.meta)
    Base.write(io, img.data)
end

read(io::IO, ::Type{RawImageHeader}) = RawImageHeader((fieldnames(RawImageHeader) .|> fieldtype $ RawImageHeader .|> MRD.unsafe_read $ io)...)

function read(io::IO, ::Type{Image})
    header = MRD.read(io, RawImageHeader)
    meta = MRD.read(io, MetaDict)
    data_type = datatype_to_type[header.data_type]

    data = Array{data_type,4}(undef, header.matrix_size[1], header.matrix_size[2], header.matrix_size[3], header.channels)
    Base.read!(io, data)
    return Image(ImageHeader(header), data, meta)

end

Base.precompile(Tuple{Type{RawImageHeader},Image})   # time: 0.089653455


text:
module EMLHeatex
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...heat_exchangers.EMLHEX_Engine
	import EMLtypes.length
	include("Heatex/Heatex_Basic.jl")
	include("Heatex/Heatex_LMTD.jl")
	include("Heatex/Heatex_NTU.jl")
end


text:
using OpenGene
using OpenGene.Reference
using ArgParse

function main(args)

    s = ArgParseSettings(description = "generate a fusion CSV file for given gene list")

    @add_arg_table s begin
        "--genes",  "-g"
            help = "the file giving gene list (separated by space, tab or line break)"
            required = true
        "--ref",  "-r"
            help = "the reference genome to use (hg19/hg37/grch38)"
            required = true
        "--fusion",   "-f"
            help = "the fusion file to output"
            required = true
    end

    options = parse_args(s) # the result is a Dict{String,Any}
    detect(options["ref"], options["genes"], options["fusion"])
end

function detect(ref, genes, fusion)
    ref = lowercase(ref)
    if !(ref in ["hg19", "grch37", "grch38"])
        println("reference <-f/--ref> should be one of hg19/grch37/grch38")
        return
    end
    input = open(genes)
    output = open(fusion, "w")
    genestr = readstring(input)
    allgenes = split(genestr)

    println("loading gencode " * ref)
    index = gencode_load(ref)
    println("loading done")

    processed_genes = []

    passed_num = 0
    for genestr in allgenes
        genelist = split(genestr, '-')
        for genename in genelist
            if genename in processed_genes
                #println("----skipped ", genename, ", already processed")
                continue
            end
            gene_attrs = split(genename, '_')
            gene = gene_attrs[1]
            transcript = ""
            if length(gene_attrs)>1
                transcript = gene_attrs[2]
            end
            gene_matches = gencode_genes(index, gene)
            if length(gene_matches)==0
                #println("----skipped ", genename, ", gene name not found in gencode")
                continue
            end
            println(genename)
            g = gene_matches[1]

            # get the correct transcript if it is specified
            target_tr = nothing
            if transcript!=""
                for tr in g.transcripts
                    if startswith(tr.id, transcript)
                        target_tr = tr
                        break
                    end
                end
                if target_tr == nothing
                    #println("----skipped ", genename, ", transcript not found")
                end
            end

            # get the 001 transcript
            if target_tr == nothing
                for tr in g.transcripts
                    tname = tr.attributes["transcript_name"]
                    if endswith(tname, "001\"")
                        target_tr = tr
                        break
                    end
                end
            end

            # get the longest transcript if it is not specified
            if target_tr == nothing
                longest = 0
                for tr in g.transcripts
                    if length(tr.exons)>longest
                        longest = length(tr.exons)
                        target_tr = tr
                    end
                end
            end

            write(output, ">", gene, "_", target_tr.id, ",", g.chr, ":", string(g.start_pos), "-", string(g.end_pos), "\n")

            exons = target_tr.exons
            id = 1
            for e in exons
                write(output, string(id), ",", string(e.start_pos), ",", string(e.end_pos), "\n")
                id += 1
            end
            passed_num += 1
            write(output, "\n")
            push!(processed_genes, genename)
        end
    end
    println("successfully processed genes: ", passed_num)
    flush(output)
end

# REPL is not supported

main(ARGS)

max_stars_repo_name: nvolkovaGEL/GeneFuse

max_issues_repo_name: nvolkovaGEL/GeneFuse

max_forks_repo_name: nvolkovaGEL/GeneFuse


text:
# Invoke with
# `julia --project=deps deps/build_local.jl [path-to-enzyme]`

# the pre-built Enzyme_jll might not be loadable on this platform
Enzyme_jll = Base.UUID("7cc45869-7501-5eee-bdea-0790c847d4ef")

using Pkg, Scratch, Preferences, Libdl

# 1. Get a scratch directory
scratch_dir = get_scratch!(Enzyme_jll, "build")
isdir(scratch_dir) && rm(scratch_dir; recursive=true)

source_dir = nothing
branch = nothing
if length(ARGS) == 2
    @assert ARGS[1] == "--branch"
    branch = ARGS[2]
    source_dir = nothing
elseif length(ARGS) == 1
    source_dir = ARGS[1]
end

if branch === nothing
    branch = "main"
end

if source_dir === nothing
    scratch_src_dir = get_scratch!(Enzyme_jll, "src")
    cd(scratch_src_dir) do
        if !isdir("Enzyme")
            run(`git clone https://github.com/wsmoses/Enzyme`)
        end
        run(`git -C Enzyme fetch`)
        run(`git -C Enzyme checkout origin/$(branch)`)
    end
    source_dir = joinpath(scratch_src_dir, "Enzyme", "enzyme")
end

# 2. Ensure that an appropriate LLVM_full_jll is installed
Pkg.activate(; temp=true)
llvm_assertions = try
    cglobal((:_ZN4llvm24DisableABIBreakingChecksE, Base.libllvm_path()), Cvoid)
    false
catch
    true
end
LLVM = if llvm_assertions
    Pkg.add(name="LLVM_full_assert_jll", version=Base.libllvm_version)
    using LLVM_full_assert_jll
    LLVM_full_assert_jll
else
    Pkg.add(name="LLVM_full_jll", version=Base.libllvm_version)
    using LLVM_full_jll
    LLVM_full_jll
end
LLVM_DIR = joinpath(LLVM.artifact_dir, "lib", "cmake", "llvm")
LLVM_VER_MAJOR = Base.libllvm_version.major

# Build!
@info "Building" source_dir scratch_dir LLVM_DIR
run(`cmake -DLLVM_DIR=$(LLVM_DIR) -DENZYME_EXTERNAL_SHARED_LIB=ON -B$(scratch_dir) -S$(source_dir)`)
run(`cmake --build $(scratch_dir) --parallel $(Sys.CPU_THREADS) -t Enzyme-$(LLVM_VER_MAJOR)`)

# Discover built libraries
built_libs = filter(readdir(joinpath(scratch_dir, "Enzyme"))) do file
    endswith(file, ".$(Libdl.dlext)") && startswith(file, "lib")
end
lib_path = joinpath(scratch_dir, "Enzyme", only(built_libs))
isfile(lib_path) || error("Could not find library $lib_path in build directory")

# Tell Enzyme_jll to load our library instead of the default artifact one
set_preferences!(
    joinpath(dirname(@__DIR__), "LocalPreferences.toml"),
    "Enzyme_jll",
    "libEnzyme_path" => lib_path;
    force=true,
)


text:
"""
`StokesSolution`

Holds the data for the solution to a Stokes problem.

### Fields

* u
* v
* p
* u_analytic
* vTrue
* pTrue
* mesh
* trueknown
* errors
* converrors

"""
type StokesSolution{T,N,U,V,P,E,C} <: AbstractNoTimeSolution{T,N}
  u::Array{T,N}
  v::Array{T,N}
  p::Array{T,N}
  u_analytic::U
  vTrue::V
  pTrue::P
  mesh::FDMMesh
  trueknown::Bool
  errors::E
  converrors::C
end

  StokesSolution(u,v,p,u_analytic,vTrue,pTrue,mesh,trueknown;errors=nothing,converrors=nothing) =
  StokesSolution(u,v,p,u_analytic,vTrue,pTrue,mesh,trueknown,errors,converrors)


text: const TISSUES = ["DLD1", "GBM", "HCT116", "HELA", "RPE1"]