text:
using BinaryBuilder, BinaryBuilder.BinaryBuilderBase, BinaryBuilder.Wizard
using GitHub, Test, VT100, Sockets, HTTP, SHA, Tar
import Pkg: PackageSpec

import BinaryBuilder.BinaryBuilderBase: available_gcc_builds, available_llvm_builds, getversion

function with_wizard_output(f::Function, state, step_func::Function)
    # Create fake terminal to communicate with BinaryBuilder over
    pty = VT100.create_pty(false)
    state.ins = Base.TTY(pty.slave)
    state.outs = Base.TTY(pty.slave)

    # Immediately start reading in off of `state.outs`
    out_buff = PipeBuffer()
    reader_task = @async begin
        while isopen(pty.master)
            z = String(readavailable(pty.master))

            # Un-comment this to figure out what on earth is going wrong
            # print(z)
            write(out_buff, z)
        end
    end

    # Start the wizard poppin' off
    wizard_task = @async begin
        try
            step_func(state)
        catch e
            bt = catch_backtrace()
            Base.display_error(stderr, e, bt)

            # If this fails, panic
            Test.@test false
        end
    end

    f(pty.master, out_buff)

    # Wait for the wizard to finish
    wait(wizard_task)

    # Once that's done, kill the reader task
    close(pty.master)
    wait(reader_task)
end

# Test the download stage
r = HTTP.Router()
io = IOBuffer()
Tar.create(joinpath(build_tests_dir, "libfoo"), pipeline(`gzip -9`, io))
libfoo_tarball_data = take!(io)
libfoo_tarball_hash = bytes2hex(sha256(libfoo_tarball_data))
function serve_tgz(req)
    HTTP.Response(200, libfoo_tarball_data)
end
HTTP.@register(r, "GET", "/*/source.tar.gz", serve_tgz)
port = -1
server = Sockets.TCPServer()
# Try to connect to different ports, in case one is busy.  Important in case we
# have multiple parallel builds.
available_ports = 14444:14544
for i in available_ports
    try
        # Update the global server to shut it down when we are done with it.
        global server = Sockets.listen(Sockets.InetAddr(Sockets.localhost, i))
    catch e
        if e isa Base.IOError
            if i == last(available_ports)
                # Oh no, this was our last attempt
                error("No more ports available for the HTTP server")
            end
            # If the port is busy, try the next one
            continue
        else
            rethrow(e)
        end
    end
    # All looks good, update the global `port` and start the server
    global port = i
    @async HTTP.serve(r, Sockets.localhost, port; server=server, verbose=false)
    break
end

function readuntil_sift(io::IO, needle)
    # N.B.: This is a terrible way to do this and works around the fact that our `IOBuffer`
    # does not block. It works fine here, but do not copy this to other places.
    needle = codeunits(needle)
    buffer = zeros(UInt8, length(needle))
    all_buffer = UInt8[]
    while isopen(io)
        new_c = read(io, 1)
        append!(all_buffer, new_c)
        if isempty(new_c)
            # We need to wait for more data, sleep for a bit
            sleep(0.01)
            continue
        end

        buffer = [buffer[2:end]; new_c]
        if !any(buffer .!= needle)
            return all_buffer
        end
    end
    return nothing
end

function call_response(ins, outs, question, answer; newline=true)
    @assert readuntil_sift(outs, question) !== nothing
    # Because we occasionally are dealing with things that do strange
    # stdin tricks like reading raw stdin buffers, we sleep here for safety.
    sleep(0.1)
    print(ins, answer)
    if newline
        println(ins)
    end
end

@testset "Wizard - Obtain source" begin
    state = Wizard.WizardState()
    # Use a non existing name
    with_wizard_output(state, Wizard.get_name_and_version) do ins, outs
        # Append "_jll" to the name and make sure this is automatically removed
        call_response(ins, outs, "Enter a name for this project", "libfoobarqux_jll")
        call_response(ins, outs, "Enter a version number", "1.2.3")
    end
    @test state.name == "libfoobarqux"
    @test state.version == v"1.2.3"
    state.name = nothing
    # Use an existing name, choose a new one afterwards
    with_wizard_output(state, Wizard.get_name_and_version) do ins, outs
        call_response(ins, outs, "Enter a name for this project", "libcurl")
        call_response(ins, outs, "Choose a new project name", "y")
        call_response(ins, outs, "Enter a name for this project", "libfoobarqux")
    end
    @test state.name == "libfoobarqux"
    @test state.version == v"1.2.3"
    state.name = nothing
    # Use an existing name, confirm the choice
    with_wizard_output(state, Wizard.get_name_and_version) do ins, outs
        call_response(ins, outs, "Enter a name for this project", "libcurl")
        call_response(ins, outs, "Choose a new project name", "N")
    end
    @test state.name == "libcurl"
    @test state.version == v"1.2.3"
end

# Set the state up
function step2_state()
    state = Wizard.WizardState()
    state.step = :step2
    state.platforms = [Platform("x86_64", "linux")]

    return state
end

@testset "Wizard - Downloading" begin
    state = step2_state()
    with_wizard_output(state, Wizard.step2) do ins, outs
        call_response(ins, outs, "Please enter a URL", "http://127.0.0.1:$(port)/a/source.tar.gz")
        call_response(ins, outs, "Would you like to download additional sources", "N")
        call_response(ins, outs, "Do you require any (binary) dependencies", "N")

        call_response(ins, outs, "Enter a name for this project", "libfoo")

        # Test bad version number detection
        call_response(ins, outs, "Enter a version number", "parse me, I dare you")
        call_response(ins, outs, "Enter a version number", "1.0.0")

        # Compiler
        call_response(ins, outs, "Do you want to customize the set of compilers?", "Y")
        call_response(ins, outs, "Select compilers for the project", "ad")
        call_response(ins, outs, "Select the preferred GCC version", "\r")
        call_response(ins, outs, "Select the preferred LLVM version", "\e[B\e[B\e[B\r")
    end
    # Check that the state is modified appropriately
    @test state.source_urls == ["http://127.0.0.1:$(port)/a/source.tar.gz"]
    @test getfield.(state.source_files, :hash) == [libfoo_tarball_hash]
    @test Set(state.compilers) == Set([:c, :rust, :go])
    @test state.preferred_gcc_version == getversion(available_gcc_builds[1])
    # The default LLVM shard is the latest one, and above we pressed three times
    # arrow down in the reverse order list.
    @test state.preferred_llvm_version == getversion(available_llvm_builds[end-3])

    # Test two tar.gz download
    state = step2_state()
    with_wizard_output(state, Wizard.step2) do ins, outs
        call_response(ins, outs, "Please enter a URL", "http://127.0.0.1:$(port)/a/source.tar.gz")
        call_response(ins, outs, "Would you like to download additional sources", "Y")
        call_response(ins, outs, "Please enter a URL", "http://127.0.0.1:$(port)/b/source.tar.gz")
        call_response(ins, outs, "Would you like to download additional sources", "N")
        call_response(ins, outs, "Do you require any (binary) dependencies", "N")

        call_response(ins, outs, "Enter a name for this project", "libfoo")
        call_response(ins, outs, "Enter a version number", "1.0.0")

        call_response(ins, outs, "Do you want to customize the set of compilers?", "N")
    end
    # Check that the state is modified appropriately
    @test state.source_urls == [
        "http://127.0.0.1:$(port)/a/source.tar.gz",
        "http://127.0.0.1:$(port)/b/source.tar.gz",
    ]
    @test getfield.(state.source_files, :hash) == [
        libfoo_tarball_hash,
        libfoo_tarball_hash,
    ]

    # Test download/install with a broken symlink that used to kill the wizard
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/183
    state = step2_state()
    with_wizard_output(state, Wizard.step2) do ins, outs
        call_response(ins, outs, "Please enter a URL", "https://github.com/staticfloat/small_bin/raw/d846f4a966883e7cc032a84acf4fa36695d05482/broken_symlink/broken_symlink.tar.gz")
        call_response(ins, outs, "Would you like to download additional sources", "N")
        call_response(ins, outs, "Do you require any (binary) dependencies", "N")

        call_response(ins, outs, "Enter a name for this project", "broken_symlink")
        call_response(ins, outs, "Enter a version number", "1.0.0")

        call_response(ins, outs, "Do you want to customize the set of compilers?", "N")
    end

    # Test failure to resolve a dependency
    state = step2_state()
    @test_logs (:warn, r"Unable to resolve iso_codez_jll") match_mode=:any with_wizard_output(state, Wizard.step2) do ins, outs
        call_response(ins, outs, "Please enter a URL", "http://127.0.0.1:$(port)/a/source.tar.gz")
        call_response(ins, outs, "Would you like to download additional sources", "N")
        call_response(ins, outs, "Do you require any (binary) dependencies", "Y")

        call_response(ins, outs, "Enter JLL package name:", "ghr_jll")
        call_response(ins, outs, "Would you like to provide additional dependencies?", "Y")
        # Test auto-JLL suffixing
        call_response(ins, outs, "Enter JLL package name:", "Zlib")
        call_response(ins, outs, "Would you like to provide additional dependencies?", "Y")

        # Test typo detection
        call_response(ins, outs, "Enter JLL package name:", "iso_codez_jll")
        call_response(ins, outs, "Unable to resolve", "N")

        call_response(ins, outs, "Enter a name for this project", "check_deps")
        call_response(ins, outs, "Enter a version number", "1.0.0")

        call_response(ins, outs, "Do you want to customize the set of compilers?", "N")
    end
    @test length(state.dependencies) == 2
    @test any([BinaryBuilder.getname(d) == "ghr_jll" for d in state.dependencies])
    @test any([BinaryBuilder.getname(d) == "Zlib_jll" for d in state.dependencies])
end

# Dump the tarball to disk so that we can use it directly in the future
tempspace = tempname()
mkdir(tempspace)
libfoo_tarball_path = joinpath(tempspace, "source.tar.gz")
open(f -> write(f, libfoo_tarball_data), libfoo_tarball_path, "w")

function step3_state()
    state = Wizard.WizardState()
    state.step = :step34
    state.platforms = [Platform("x86_64", "linux")]
    state.source_urls = ["http://127.0.0.1:$(port)/a/source.tar.gz"]
    state.source_files = [BinaryBuilder.SetupSource{ArchiveSource}(libfoo_tarball_path, libfoo_tarball_hash, "")]
    state.name = "libfoo"
    state.version = v"1.0.0"
    state.dependencies = Dependency[]
    state.compilers = [:c]
    state.preferred_gcc_version = getversion(available_gcc_builds[1])
    state.preferred_llvm_version = getversion(available_llvm_builds[end])

    return state
end

function step3_test(state)
    @test length(state.files) == 2
    @test "lib/libfoo.so" in state.files
    @test "bin/fooifier" in state.files

    libfoo_idx = findfirst(state.files .== "lib/libfoo.so")
    fooifier_idx = findfirst(state.files .== "bin/fooifier")
    @test state.file_kinds[libfoo_idx] == :library
    @test state.file_kinds[fooifier_idx] == :executable
    @test state.file_varnames[libfoo_idx] == :libfoo
    @test state.file_varnames[fooifier_idx] == :fooifier
end

@testset "Wizard - Building" begin
    function succcess_path_call_response(ins, outs)
        output = readuntil_sift(outs, "Build complete")
        if contains(String(output), "Warning:")
            close(ins)
            return false
        end
        call_response(ins, outs, "Would you like to edit this script now?", "N")
        call_response(ins, outs, "d=done, a=all", "ad"; newline=false)
        call_response(ins, outs, "lib/libfoo.so", "libfoo")
        call_response(ins, outs, "bin/fooifier", "fooifier")
        return true
    end

    # Test step3 success path
    state = step3_state()
    with_wizard_output(state, Wizard.step34) do ins, outs
        call_response(ins, outs, "\${WORKSPACE}/srcdir", """
        cd libfoo
        make install
        exit
        """)
        @test succcess_path_call_response(ins, outs)
    end
    @test state.history == """
    cd \$WORKSPACE/srcdir
    cd libfoo
    make install
    exit
    """
    step3_test(state)

    # Step 3 failure path (no binary in destdir -> return to build)
    state = step3_state()
    with_wizard_output(state, Wizard.step34) do ins, outs
        # Don't build anything
        call_response(ins, outs, "\${WORKSPACE}/srcdir", "exit")
        call_response(ins, outs, "Would you like to edit this script now?", "N")

        # Return to build environment
        call_response(ins, outs, "Return to build environment", "\r", newline=false)
        call_response(ins, outs, "\${WORKSPACE}/srcdir", """
        cd libfoo
        make install
        exit
        """)

        @test succcess_path_call_response(ins, outs)
    end
    @test state.history == """
    cd \$WORKSPACE/srcdir
    exit
    cd \$WORKSPACE/srcdir
    cd libfoo
    make install
    exit
    """
    step3_test(state)

    # Step 3 failure path (no binary in destdir -> retry with a clean build environment)
    state = step3_state()
    with_wizard_output(state, Wizard.step34) do ins, outs
        # Don't build anything
        call_response(ins, outs, "\${WORKSPACE}/srcdir", "exit")
        call_response(ins, outs, "Would you like to edit this script now?", "N")

        # Clean environment
        call_response(ins, outs, "Return to build environment", "\e[B\r")
    end
    @test state.step == :step3

    # Step 3 with a failing script
    state = step3_state()
    with_wizard_output(state, Wizard.step34) do ins, outs
        # Build ok, but then indicate a failure
        call_response(ins, outs, "\${WORKSPACE}/srcdir", """
        cd libfoo
        make install
        exit 1
        """)

        @test readuntil_sift(outs, "Warning:") !== nothing

        @test succcess_path_call_response(ins, outs)
    end
    step3_test(state)

    # Step 3 dependency download
    state = step3_state()
    state.dependencies = [Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))]
    with_wizard_output(state, Wizard.step34) do ins, outs
        call_response(ins, outs, "\${WORKSPACE}/srcdir", """
        if [[ ! -f \${libdir}/libz.\${dlext} ]]; then
            echo "ERROR: Could not find libz.\${dlext}" >&2
            exit 1
        fi
        cd libfoo
        make install
        exit
        """)
        @test succcess_path_call_response(ins, outs)
    end

    # Step 3 - `bb add`
    state = step3_state()
    state.dependencies = [Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))]
    with_wizard_output(state, Wizard.step34) do ins, outs
        call_response(ins, outs, "\${WORKSPACE}/srcdir", """
        if [[ ! -f \${libdir}/libz.\${dlext} ]]; then
            echo "ERROR: Could not find libz.\${dlext}" >&2
            exit 1
        fi
        bb add Xorg_xorgproto_jll
        if [[ ! -d \${includedir}/X11 ]]; then
            echo "ERROR: Could not find include/X11" >&2
            exit 1
        fi
        bb add Zlib_jll
        cd libfoo
        make install
        exit
        """)
        @test succcess_path_call_response(ins, outs)
    end
end

function step5_state(script)
    state = step3_state()
    state.history = script
    state.files = ["lib/libfoo.so","bin/fooifier"]
    state.file_kinds = [:library, :executable]
    state.file_varnames = [:libfoo, :fooifier]
    state
end

@testset "Wizard - Generalizing" begin
    # Check that with a failing script, step 5 rejects,
    # even if all artifacts are present.
    state = step5_state("""
        cd libfoo
        make install
        exit 1
    """)

    with_wizard_output(state, state->Wizard.step5_internal(state, first(state.platforms))) do ins, outs
        call_response(ins, outs, "Press Enter to continue...", "\n")
        call_response(ins, outs, "How would you like to proceed?", "\e[B\e[B\r")
    end
    @test isempty(state.platforms)
end

function step7_state()
    state = step5_state("""
    cd libfoo
    make install
    exit 1
    """)
    state.patches = [PatchSource("foo.patch", "this is a patch")]
    return state
end

@testset "Wizard - Deployment" begin
    state = step7_state()
    # First, test local deployment
    mktempdir() do out_dir
        with_wizard_output(state, state->Wizard._deploy(state)) do ins, outs
            call_response(ins, outs, "How should we deploy this build recipe?", "\e[B\r")
            call_response(ins, outs, "Enter directory to write build_tarballs.jl to:", "$(out_dir)\r")
        end
        @test isfile(joinpath(out_dir, "build_tarballs.jl"))
        @test isfile(joinpath(out_dir, "bundled", "patches", "foo.patch"))
    end

    # Next, test writing out to stdout
    state = step7_state()
    with_wizard_output(state, state->Wizard._deploy(state)) do ins, outs
        call_response(ins, outs, "How should we deploy this build recipe?", "\e[B\e[B\r")
        @test readuntil_sift(outs, "Your generated build_tarballs.jl:") !== nothing
        @test readuntil_sift(outs, "name = \"libfoo\"") !== nothing
        @test readuntil_sift(outs, "make install") !== nothing
        @test readuntil_sift(outs, "LibraryProduct(\"libfoo\", :libfoo)") !== nothing
        @test readuntil_sift(outs, "ExecutableProduct(\"fooifier\", :fooifier)") !== nothing
        @test readuntil_sift(outs, "dependencies = Dependency[") !== nothing
    end
end

@testset "Wizard - state serialization" begin
    for state_generator in (Wizard.WizardState, step2_state, step3_state, step7_state)
        mktempdir() do dir
            state = state_generator()

            Wizard.save_wizard_state(state, dir)
            @test Wizard.load_wizard_state(dir; as_is=true) == state
        end
    end
end

close(server)

@testset "GitHub - authentication" begin
    withenv("GITHUB_TOKEN" => "") do
        @test Wizard.github_auth(allow_anonymous=true) isa GitHub.AnonymousAuth
    end
end


text:
function setup_model(wordembeds::Vector{Var}, charembeds::Vector{Var}, ntags::Int)
    T = eltype(wordembeds[1])
    w = Node()
    hw = lookup(wordembeds, w)

    c = Node()
    hc = lookup(charembeds, c)
    d = size(charembeds[1], 1)
    hc = Conv1D(T,5,d,5d,2,1)(hc)
    hc = max(hc, 2)
    hc = resize(hc, batchsize(w))

    h = concat(1, hw, hc)
    d = size(wordembeds[1],1) + 5size(charembeds[1],1)
    dh = 300
    h = Conv1D(T,5,d,dh,2,1)(h)
    h = relu(h)

    #x = Node(dropout, x, 0.3)
    #x = Node(Conv1D(T,5,dh,dh,2,1), x)
    #x = Node(relu, x)
    #x = Node(Linear(T,dh,ntags), x)
    Graph(input=(w,c), output=h)
end

#=
function setup_posembeds{T}(::Type{T}, dim::Int, len::Int)
    embeds = zeros(T, dim, len)
    for p = 1:len
        for i = 1:2:dim
            x = p / 10000^((i-1)/dim)
            embeds[i,p] = sin(x)
            embeds[i+1,p] = cos(x)
        end
    end
    embeds
end
=#

ext: jl

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using OrdinaryDiffEq, Test, DiffEqBase, DiffEqCallbacks, RecursiveArrayTools

u0 = ones(2,2)
f = function (du,u,p,t)
  du[1,:] = u[2,:]
  du[2,:] = -u[1,:]
end
prob = ODEProblem(f,u0,(0.0,100.0))

function g(resid,u,p)
  resid[1] = u[2]^2 + u[1]^2 - 2
  resid[2] = u[3]^2 + u[4]^2 - 2
  resid[3] = 0
  resid[4] = 0
end

g_t(resid,u,p,t) = g(resid,u,p)

isautonomous(p::ManifoldProjection{autonomous,NL}) where {autonomous,NL} = autonomous

sol = solve(prob,Vern7())
@test !(sol[end][1]^2 + sol[end][2]^2 ≈ 2)

# autodiff=true
cb = ManifoldProjection(g)
@test isautonomous(cb.affect!)
solve(prob,Vern7(),callback=cb)
@time sol=solve(prob,Vern7(),callback=cb)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2

cb_t = ManifoldProjection(g_t)
@test !isautonomous(cb_t.affect!)
solve(prob,Vern7(),callback=cb_t)
@time sol_t = solve(prob,Vern7(),callback=cb_t)
@test sol_t.u == sol.u && sol_t.t == sol.t

# autodiff=false
cb_false = ManifoldProjection(g, nlsolve=DiffEqCallbacks.NLSOLVEJL_SETUP(autodiff=false))
@test isautonomous(cb_false.affect!)
solve(prob,Vern7(),callback=cb_false)
sol=solve(prob,Vern7(),callback=cb_false)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2

cb_t_false = ManifoldProjection(g_t, nlsolve=DiffEqCallbacks.NLSOLVEJL_SETUP(autodiff=false))
@test !isautonomous(cb_t_false.affect!)
solve(prob,Vern7(),callback=cb_t_false)
sol_t = solve(prob,Vern7(),callback=cb_t_false)
@test sol_t.u == sol.u && sol_t.t == sol.t

# test array partitions
u₀ = ArrayPartition(ones(2), ones(2))
prob = ODEProblem(f, u₀, (0.0, 100.0))

sol = solve(prob,Vern7(),callback=cb)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2

sol = solve(prob,Vern7(),callback=cb_t)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2

# does not work since Calculus.jl (on which NLsolve.jl depends)
# implements only Jacobians of vectors
sol = solve(prob,Vern7(),callback=cb_false)
sol[end][1]^2 + sol[end][2]^2 ≈ 2

sol = solve(prob,Vern7(),callback=cb_t_false)
sol[end][1]^2 + sol[end][2]^2 ≈ 2


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m5_opt_rho1_68(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m5_opt_rho1_68!(A_copy)
end

@inline function exp_ps_m5_opt_rho1_68!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=0.3689088315142506
    coeff2=0.993511502899774
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.3689088315142506
    coeff2=0.993511502899774
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.06095437116133005
    coeff2=0.35179926100509473
    coeff3=1.2664069151339148
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.5306796402178436
    coeff2=1.2516399114571966
    coeff3=0.08935374652915735
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.014130431015307159
    coeff2=-0.030997495252697738
    coeff3=0.33768096635087375
    coeff4=0.9844379024704415
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=1.0138346120944015
    coeff2=0.3055672291311765
    coeff3=0.045499149053109256
    coeff4=0.0014555164226796927
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-1.7579891447774443e-7
    coeff2=0.0013217237148109482
    coeff3=0.0061860694149101746
    coeff4=-0.03360537533961365
    coeff5=1.0393506245661515
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=5.01359878707017e-6
    coeff2=-1.3147702734172197e-7
    coeff3=2.1793867756624317e-8
    coeff4=4.212389951082164e-10
    coeff5=1.9503051245595482e-11
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.002218807643492532
    coeff2=-0.02673358365584956
    coeff3=-0.13331176967937897
    coeff4=0.9627486812176257
    coeff5=0.7011206850897305
    coeff6=4.940845221165412e-5
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.05934666479458639
    coeff2=0.003857202961892145
    coeff3=0.0005189095187953901
    coeff4=3.843013696548252e-5
    coeff5=2.1244028410464186e-6
    coeff6=1.0000000004106206
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.9526549410501626
    coeff2=0.7066204202720486
    coeff3=0.31964768881701877
    coeff4=-0.006731411055840937
    coeff5=-0.020134332996101654
    coeff6=-1.4225442860061685e-5
    coeff7=0.9991129755494675
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
Sorry_OverdubInstead_Is_Deprecated() = error("Usage of `OverdubInstead` should be replaced by `Cassette.recurse`; see `Cassette.recurse` for details.")

Base.@deprecate_binding(OverdubInstead, Sorry_OverdubInstead_Is_Deprecated, false,
string("; the `overdub` <--> `execute` cycle has been removed from Cassette, and",
       " usage of `OverdubInstead` has been replaced with `Cassette.recurse`;",
       " see `Cassette.recurse` for details."))

struct Sorry_Execute_Is_Deprecated end

Base.@deprecate_binding(execute, Sorry_Execute_Is_Deprecated, false,
string("; the `overdub` <--> `execute` cycle has been removed from Cassette. In",
       " most cases, `execute` calls can now be replaced with `overdub`; in",
       " general, however, users may need to do some extra refactoring to",
       " properly migrate to the new overdubbing model. Feel free to open a",
       " Cassette issue for help with migration."))

Base.@deprecate(canoverdub(args...), canrecurse(args...))


text:
sum_bcast(x, y) = sum(x .+ y)

@testset "special bcast" begin
    for args in [
        (rand(3, 4), rand(3)),
        (rand(3, 4), rand(3, 1)),
        (rand(3, 4), rand(1, 4)),
        (rand(3), rand(3, 4)),
        (rand(3, 1), rand(3, 4)),
        (rand(1, 4), rand(3, 4)),
    ]
        val, g = grad(sum_bcast, args...; static=false)
        for i=1:length(args)
            @test size(g[i]) == size(args[i])
        end
    end
end

mutable struct Linear{T}
    W::AbstractArray{T,2}
    b::AbstractArray{T}
end

forward(m::Linear, X) = m.W * X

loss(m::Linear, X) = sum(forward(m, X))

@testset "grad: structs" begin

    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]
    val, g = grad(loss, args...)

    @test val == loss(args...)

    play!(g.tape)
    val, g = grad(loss, args...)
    @test val == loss(args...)

end

HESS = randn(3,3)
hessian_fun(x) = x'*(HESS*x)
hessian_fun2(x) = 0.5*x'*(HESS*x)

@testset "grad: adjoint" begin
    H = randn(3,3)
    x = rand(3)
    val, g = grad(hessian_fun2, x)
    @test val isa Real
end

@testset "grad: compiled" begin
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]

    # static vs. dynamic
    val1, g1 = grad(loss, args...; static=false)
    grad(loss, args...; static=true)
    val2, g2 = grad(loss, args...; static=true)
    @test val1 == val2
    @test g1[1][(:W,)] == g2[1][(:W,)]

    # compiled vs. non-compiled
    tape = g1.tape
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]

    play!(tape, args...; use_compiled=true)
    last_val1 = getvalue(tape[end])
    play!(tape, args...; use_compiled=false)
    last_val2 = getvalue(tape[end])
    @test last_val1 == last_val2

    # (* -> mul!) for mixed array-scalar vars
    x = rand(3)
    val1, g1 = grad(hessian_fun, x)  # interpreted
    val2, g2 = grad(hessian_fun, x)  # compiled
    @test val1 == val2
    @test g1[1] == g2[1]
end


text:
using Documenter, ApproxInferenceProblems

makedocs(
    modules = [ApproxInferenceProblems],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Johanni Brea",
    sitename = "ApproxInferenceProblems.jl",
    pages = Any["index.md"],
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(repo = "github.com/jbrea/ApproxInferenceProblems.jl.git", push_preview = true)


text:
function check_iter(
    current_iter::Int64;
    max_iter::Int64 = 100,
    verbosity::Int64 = 2
    )
    if verbosity > 1
        println("Iteration: #", current_iter)
    end
    if current_iter >= max_iter
    println(
            "Maximum number of iterations reached after ",
            Int(current_iter),
            " iterations.",
        )
        return true
    else
        return false
    end
end

function check_time(
    start_time::Float64;
    max_time::Int64 = 1000,
    verbosity::Int64 = 2
    )
    current_time = time()-start_time
    if current_time>= max_time
    println(
            "Max time reached after ",
            current_time,
            " seconds.",
        )
        return true
    else
        return false
    end
end

function check_f_tol_rel!(
    examinees::Vector{<:AbstractExaminee},
    old_likelihood::Float64;
    f_tol_rel::Float64 = 0.000001,
    verbosity::Int64 = 2
    )
    new_likelihood = likelihood(examinees)
    f_rel = abs((new_likelihood - old_likelihood) / old_likelihood)
    if verbosity > 1
        println("Likelihood: ", new_likelihood)
    end
    if f_rel < f_tol_rel
        println(
            "Likelihood maximum relative tolerance reached."
        )
        println("Likelihood: ", new_likelihood)
        old_likelihood = copy(new_likelihood)
        return true
    else
        old_likelihood = copy(new_likelihood)
        return false
    end
end

function check_f_tol_rel!(
    latents::Vector{<:AbstractLatent},
    old_likelihood::Float64;
    f_tol_rel::Float64 = 0.000001,
    verbosity::Int64 = 2
    )
    new_likelihood = _likelihood(latents)
    f_rel = abs((new_likelihood - old_likelihood) / old_likelihood)
    if verbosity > 1
        println("Likelihood: ", new_likelihood)
    end
    if f_rel < f_tol_rel
        println(
            "Likelihood maximum relative tolerance reached."
        )
        println("Likelihood: ", new_likelihood)
        old_likelihood = copy(new_likelihood)
        return true
    else
        old_likelihood = copy(new_likelihood)
        return false
    end
end

function check_f_tol_rel!(
    new_likelihood::Float64,
    old_likelihood::Float64;
    f_tol_rel::Float64 = 0.01,
    verbosity::Int64 = 2
    )
    f_rel = abs((new_likelihood - old_likelihood) / old_likelihood)
    if verbosity > 1
        println("Likelihood: ", new_likelihood)
    end
    if f_rel < f_tol_rel
        println(
            "Likelihood maximum relative tolerance reached."
            )
        println("Likelihood: ", new_likelihood)
        old_likelihood = copy(new_likelihood)
        return true
    else
        old_likelihood = copy(new_likelihood)
        return false
    end
end

function check_x_tol_rel!(
    new_pars::Matrix{Float64},
    old_pars::Matrix{Float64};
    x_tol_rel::Float64 = 0.0001,
    verbosity::Int64 = 2
)
    delta_pars =  maximum(abs.((new_pars - old_pars) ./ old_pars))
    if verbosity > 1
        println("x_rel max: ", delta_pars)
    end
    if delta_pars <= x_tol_rel
        println(
            "Item parameters maximum relative tolerance reached."
            )
        println("x_rel max: ", delta_pars)
        old_pars .= new_pars
        return true
    else
        old_pars .= new_pars
        return false
    end
end

function check_x_tol_rel!(
    parameters::Vector{<:AbstractParameters},
    old_pars::Matrix{Float64};
    x_tol_rel::Float64 = 0.0001,
    verbosity::Int64 = 2
)
    new_pars = hcat(_get_parameters_vals.(parameters)...)
    delta_pars =  maximum(abs.((new_pars - old_pars) ./ old_pars))
    if verbosity > 1
        println("x_rel max: ", delta_pars)
    end
    if delta_pars <= x_tol_rel
        println(
            "Item parameters maximum relative tolerance reached."
        )
        println("x_rel max: ", delta_pars)
        old_pars .= new_pars
        return true
    else
        old_pars .= new_pars
        return false
    end
end

function check_x_tol_rel!(
    items::Vector{<:AbstractItem},
    old_pars::Matrix{Float64};
    x_tol_rel::Float64 = 0.0001,
    verbosity::Int64 = 2
)
    new_pars = get_parameters_vals(items)
    delta_pars =  maximum(abs.((new_pars - old_pars) ./ old_pars))
    if verbosity > 1
        println("x_rel max: ", delta_pars)
    end
    if delta_pars <= x_tol_rel
        println(
            "Item parameters maximum relative tolerance reached."
            )
        old_pars .= new_pars
        return true
    else
        old_pars .= new_pars
        return false
    end
end

function check_x_tol_rel!(
    latents::Vector{<:AbstractLatent},
    old_latents::Matrix{Float64};
    x_tol_rel::Float64 = 0.0001,
    verbosity::Int64 = 2
)
    new_latents = hcat(_get_latents.(latents)...)
    delta_latents =  maximum(abs.((new_latents - old_latents) ./ old_latents))
    if verbosity > 1
        println("x_rel max: ", delta_pars)
    end
    if delta_latents <= x_tol_rel
        println(
            "Item parameters maximum relative tolerance reached."
            )
        println("x_rel max: ", delta_pars)
        old_latents .= new_latents
        return true
    else
        old_latents .= new_latents
        return false
    end
end


text:
using .Gadfly

using NeXLMatrixCorrection
using Colors

"""
    Gadfly.plot(
        dbkrs::AbstractArray{DBKRatio};
        mc::Type{<:MatrixCorrection} = XPP,
        fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
        palette = NeXLPalette,
        style::Symbol = :Ratio | :XY
    )
"""
function Gadfly.plot(
    dbkrs::AbstractArray{DBKRatio};
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
    style::Symbol = :Ratio
)
    if style==:Ratio
        plot_ratio(dbkrs, mc, fc, cc)
    else
        plot_xy(dbkrs, mc, fc)
    end
end

function plot_ratio(
    dbkrs::AbstractArray{DBKRatio},
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
)
  df = asa(DataFrame, dbkrs, withComputedKs=true, mc=mc, fc=fc, cc=cc)
  df[:, "Element"] = symbol.(element.(first.(df[:,"Lines"])))
  df[:, "Z"] = z.(element.(first.(df[:,"Lines"])))
  df[:, "Family"] = name.(df[:,"Lines"], true)
  df[:, "RU"] = map(r->uv(r["K"],r["ΔK"])/r["Kxpp"], eachrow(df))
  df[:, "Rmin"] = map(r-> value(r["RU"])-σ(r["RU"]), eachrow(df))
  df[:, "Rmax"] = map(r-> value(r["RU"])+σ(r["RU"]), eachrow(df))
  df[:,"Measurement"] = map(r->"$(name(r["Cmeas"])) at $(r["E0meas"]/1000.0) keV", eachrow(df))
  df[:,"Rand"] = rand(nrow(df)) # To plot points in a randomized order.
  filter!(r->element(first(r[:Lines])) in keys(r[:Cmeas]), df)
  sort!(df, [:Family, :E0meas, :Z, :Rand])
  # yin, yax = something(ymin, 0.9*minimum(skipmissing(df[:,"Rmin"]))), something(ymax, 1.1*maximum(skipmissing(df[:,"Rmax"])))
  plot(df, xgroup="Family", x="Element", ymin="Rmin", ymax="Rmax", y="Ratio", color="Measurement",
     # Coord.cartesian(ymin=yin, ymax=yax), Stat.x_jitter(range=1.0), # Neither seems to work...
     Geom.subplot_grid(Geom.errorbar), Stat.x_jitter(range=1000.0))
end

function plot_xy(
    dbkrs::AbstractArray{DBKRatio},
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
)
    x, y, dy, color = Float64[], Float64[], Float64[], String[]
    next = 0
    for dbkr in dbkrs
        kr, unkComp = asa(KRatio, dbkr), get(dbkr.measured, :Composition, missing)
        if hasminrequired(mc, kr.unkProps) &&
           hasminrequired(fc, kr.unkProps) && #
           hasminrequired(mc, kr.stdProps) &&
           hasminrequired(fc, kr.stdProps) && #
           (!isnothing(unkComp)) &&
           (value(unkComp[kr.element]) > 0.0) &&
           (value(kr.standard[kr.element]) > 0.0)
            try
                # Compute the k-ratio
                kc = gZAFc(kr, unkComp, mc=mc, fc=fc, cc=cc) * (value(unkComp[kr.element]) / value(kr.standard[kr.element]))
                push!(y, value(kr.kratio))
                push!(dy, σ(kr.kratio))
                push!(x, kc)
                push!(color, "$(symbol(kr.element)) $(name(shell(brightest(kr.lines)))) in $(name(unkComp)) $(kr.unkProps[:BeamEnergy]/1000.0) keV")
            catch c
                @info "Failed on $dbkr - $c"
            end
        end
    end
    abline = Geom.abline(color="red", style=:dash)
    plot(
        x = x,
        y = y,
        ymin = y .- dy,
        ymax = y .+ dy,
        color = color,
        Geom.errorbar,
        Guide.xlabel("k[$(repr(nameof(mc))[2:end])]"),
        Guide.ylabel("k[Measured]"),
        Guide.colorkey(title="Measurement"),
        abline, intercept=[0.0],slope=[1.0],
    )
end

function plot2(dbkrs::AbstractArray{DBKRatio}; palette = NeXLPalette)
    kok, dkok, color = Float64[], Float64[], Color[]
    ygroup = CharXRay[]
    matcolors = Dict{String,RGB{Float64}}()
    for dbkr in dbkrs
        kr, unkComp = asa(KRatio, dbkr), get(dbkr.measured, :Composition, missing)
        push!(kok, value(kr.kratio))
        push!(dkok, σ(kr.kratio))
        push!(ygroup, brightest(kr.lines))
        matname = "$(kr.element) in $(name(unkComp)) $(kr.unkProps[:BeamEnergy]/1000.0) keV"
        if !haskey(matcolors, matname)
            matcolors[matname] = palette[length(matcolors)+1]
        end
        push!(color, matcolors[matname])
    end
    plot(
        x = eachindex(kok),
        y = kok,
        ymin = kok .- dkok,
        ymax = kok .+ dkok,
        color = color,
        ygroup = ygroup, # Geom.errorbar, Stat.x_jitter(range=0.8),
        Geom.subplot_grid(Geom.errorbar, free_y_axis = true),
        Scale.ygroup(labels = cxr -> repr(cxr), levels = allcxrs),
        Guide.manual_color_key("Material", collect(keys(matcolors)), collect(values(matcolors))), # Guide.yrug,
        Guide.xlabel("Index"),
        Guide.ylabel("k[Measured]"),
    )
end

function plot3(
    krs::AbstractArray{DBKRatio};
    known::Union{Material,Missing} = missing,
    label::AbstractString = "Material",
    palette = NeXLPalette,
)
    mats = unique(collect(map(dbkr->get(dbkr.measured,:Composition,missing), krs)))
    allelms = sort(convert(Vector{Element}, collect(union(map(keys, mats)...))))
    elmcol = Dict(elm => palette[i] for (i, elm) in enumerate(allelms))
    xs, ymin, ymax, ygroups, colors = String[], Float64[], Float64[], Element[], Color[]
    for mat in mats
        append!(xs, [name(mat) for elm in keys(mat)])
        append!(ymin, [value(mat[elm]) - σ(mat[elm]) for elm in keys(mat)])
        append!(ymax, [value(mat[elm]) + σ(mat[elm]) for elm in keys(mat)])
        append!(colors, [elmcol[elm] for elm in keys(mat)])
        append!(ygroups, collect(keys(mat)))
    end
    plot(
        x = xs,
        ymin = ymin,
        ymax = ymax,
        color = colors,
        ygroup = ygroups,
        Geom.subplot_grid(Geom.errorbar, free_y_axis = true),
        Scale.ygroup(labels = elm -> symbol(elm), levels = allelms),
        Guide.manual_color_key("Material", symbol.(collect(keys(elmcol))), collect(values(elmcol))),
        Guide.xlabel(label),
        Guide.ylabel("Mass Fraction by Element"),
    )
end


text:
module PowerPC64
  using ..Registers
  abstract RegisterSet <: Registers.RegisterSet
  using MachO
  import Base: copy

  immutable PowerPC64Arch <: Registers.Architecture
  end
  Registers.intptr(::PowerPC64Arch) = UInt64

  # See http://refspecs.linuxfoundation.org/ELF/ppc64/PPC-elf64abi.html
  const dwarf_numbering = Dict{Int, Symbol}(
    (i => Symbol("r$i") for i = 0:31)...,
    (32+i => Symbol("f$i") for i = 0:31)...,
    64 => :cr, 65 => :fpscr, 101 => :xer, 108 => :lr,
    109 => :ctr, 110 => :pc
  )
  const inverse_dwarf = map(p->p[2]=>p[1], dwarf_numbering)
  const basic_regs = sort(collect(keys(dwarf_numbering)))
  const extended_registers = 17:32

  const gdb_numbering = Dict{Int, Symbol}(
    (i => Symbol("r$i") for i = 0:31)...,
    (32+i => Symbol("f$i") for i = 0:31)...,
    64 => :pc, 65 => :msr, 66 => :cr, 67 => :lr,
    68 => :ctr, 69 => :xer, 70 => :fpscr
  )
  const inverse_gdb = map(p->p[2]=>p[1], gdb_numbering)

  const kernel_order = [
    [Symbol("r$r") for r in 0:31]...,
    :pc, :msr, :orig_gpr, :ctr, :lr, :xer, :cr, :softe,
    :trap, :dar, :dsisr, :result
  ]
  const kernel_numbering = Dict(sym=>idx for (idx,sym) in enumerate(kernel_order))

  # Basic Register Set
  const RegT = RegisterValue{UInt64}
  const ExtendedRegT = RegisterValue{UInt128}
  @eval type BasicRegs <: RegisterSet
      $(Expr(:block, (
        :($(dwarf_numbering[i]) :: RegT) for i in basic_regs
      )...))
      BasicRegs() = new()
  end
  Registers.getarch(::BasicRegs) = PowerPC64Arch()

  import ..Registers: ip, set_ip!, set_sp!, invalidate_regs!

  ip(regs::BasicRegs) = regs.pc
  set_ip!(regs::BasicRegs, ip) = regs.pc = RegisterValue{UInt64}(ip)
  set_sp!(regs::BasicRegs, sp) = regs.r1 = RegisterValue{UInt64}(sp)
  function invalidate_regs!(regs::BasicRegs)
      for fieldi = 1:nfields(regs)
          setfield!(regs, fieldi, Registers.invalidated(getfield(regs, fieldi)))
      end
  end

  function copy(regs::BasicRegs)
      ret = BasicRegs()
      for i = 1:nfields(regs)
          setfield!(ret, i, getfield(regs, i))
      end
      ret
  end

  function Base.show(io::IO, regs::BasicRegs)
      for (i,reg) in enumerate(fieldnames(typeof(regs)))
          println(io," "^(5-length(string(reg))),reg," ",getfield(regs, i))
      end
  end

  function Registers.set_dwarf!(regs::BasicRegs, reg::Integer, value)
      (reg <= last(basic_regs)) && setfield!(regs, reg+1, RegisterValue{UInt64}(value))
  end

  function Registers.get_dwarf(regs::BasicRegs, reg::Integer)
      (reg <= last(basic_regs)) ? getfield(regs, reg+1) :
        RegisterValue{UInt64}(0, 0)
  end

  function Registers.get_dwarf(::PowerPC64Arch, regs, sym)
    get_dwarf(regs, inverse_dwarf[sym])
  end

  function Registers.set_dwarf!(::PowerPC64Arch, RC, sym, val)
    set_dwarf!(RC, inverse_dwarf[sym], val)
  end

  function Registers.get_syscallarg(::PowerPC64Arch, regs, idx)
    @assert 1 <= idx <= 6
    get_dwarf(regs, [:rdi, :rsi, :rdx, :r10, :r8, :r9][idx])
  end

end


text:
using LinearAlgebra, SparseArrays, SuiteSparse

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

A = Float64[10 1 -5;-20 3 20;5 3 5]
b = Float64[1; 2; 6]

(Al, Au, Ap) = lu(A)

@show y = Al \ b[Ap]
println()

@show x = Au \ y
println()

@show As = sparse(A)
println()

@show As \ y
println()

@show F = lu(As)
println()

@show ys = F.L \ (F.Rs .* b)
println()

@show x = F.U \ ys
println()

@show F.L * F.U == (F.Rs .* A)[F.p, F.q]
println()

@show F.L * F.U
println()

@show sparse((F.Rs .* A)[F.p, F.q])
println()

end


text:
struct ADMMOptions
  ρ::Float64
  α::Float64
  maxiter::Int64
  abstol::Float64
  reltol::Float64
end

ADMMOptions(;ρ::Float64=1.,
           α::Float64=1.,
           maxiter::Int64=1000,
           abstol::Float64=1e-4,
           reltol::Float64=1e-2) = ADMMOptions(ρ, α, maxiter, abstol, reltol)

#########################################################
# minimize  trace(S*X) - log det X + lambda*||X||_1
#########################################################

function covsel!(
  X::StridedMatrix{T},
  Z::StridedMatrix{T},
  U::StridedMatrix{T},
  Σ::StridedMatrix{T},
  λ::T;
  options::ADMMOptions = ADMMOptions(),
  penalize_diag::Bool=true
  ) where {T<:AbstractFloat}

  maxiter = options.maxiter
  ρ = options.ρ
  α = options.α
  abstol = options.abstol
  reltol = options.reltol
  γ = one(T) / ρ
  λρ = λ * γ

  p = size(Σ, 1)
  g = ProxGaussLikelihood(Symmetric(Σ))

  tmpStorage = zeros(T, (p, p))
  Zold = copy(Z)

  for iter=1:maxiter
    # x-update
    @. tmpStorage = Z - U
    prox!(g, X, tmpStorage, γ)

    # z-update with relaxation
    copyto!(Zold, Z)
    @. tmpStorage = α*X + (one(T)-α)*Z + U
    @inbounds for c=1:p
      for r=1:c-1
        t = shrink(tmpStorage[r, c], λρ)
        Z[c, r] = t
        Z[r, c] = t
      end
      Z[c, c] = penalize_diag ? shrink(tmpStorage[c, c], λρ) : tmpStorage[c, c]
    end

    # u-update
    @. U = tmpStorage - Z

    # check convergence
    r_norm = norm_diff(X, Z, 2)
    s_norm = norm_diff(Z, Zold, 2) * sqrt(ρ)
    eps_pri = p*abstol + reltol * max( norm(X), norm(Z) )
    eps_dual = p*abstol + reltol * ρ * norm(U)
    if r_norm < eps_pri && s_norm < eps_dual
      break
    end
  end
  Z
end

function covselpath(S::StridedMatrix{T},
                    λarr;
                    options::ADMMOptions = ADMMOptions(),
                    penalize_diag::Bool=true,
                    verbose::Bool=false) where {T<:AbstractFloat}
  p = size(S, 1)

  solutionpath = Array(Array{Float64, 2}, length(λarr))
  X = zeros(p, p)
  Z = zeros(p, p)
  U = zeros(p, p)

  for i=1:length(λarr)
    if verbose
      @printf("lambda = %d/%d\n", i, length(λarr))
    end
    covsel!(X, Z, U, S, λarr[i]; options=options, penalize_diag=penalize_diag)
    solutionpath[i] = copy(Z)
  end
  solutionpath
end


text:
#=
Generate a MarkovChain randomly.

@author : Daisuke Oyama

=#
import StatsBase: sample
import QuantEcon: MarkovChain

# random_markov_chain

"""
Return a randomly sampled MarkovChain instance with n states.

##### Arguments

- `n::Integer` : Number of states.

##### Returns

- `mc::MarkovChain` : MarkovChain instance.

##### Examples

```julia
julia> using QuantEcon

julia> mc = random_markov_chain(3)
Discrete Markov Chain
stochastic matrix:
3x3 Array{Float64,2}:
 0.281188  0.61799   0.100822
 0.144461  0.848179  0.0073594
 0.360115  0.323973  0.315912

```

"""
function random_markov_chain(n::Integer)
    p = random_stochastic_matrix(n)
    mc = MarkovChain(p)
    return mc
end

"""
Return a randomly sampled MarkovChain instance with n states, where each state
has k states with positive transition probability.

##### Arguments

- `n::Integer` : Number of states.

##### Returns

- `mc::MarkovChain` : MarkovChain instance.

##### Examples

```julia
julia> using QuantEcon

julia> mc = random_markov_chain(3, 2)
Discrete Markov Chain
stochastic matrix:
3x3 Array{Float64,2}:
 0.369124  0.0       0.630876
 0.519035  0.480965  0.0
 0.0       0.744614  0.255386

```

"""
function random_markov_chain(n::Integer, k::Integer)
    p = random_stochastic_matrix(n, k)
    mc = MarkovChain(p)
    return mc
end

# random_stochastic_matrix

"""
Return a randomly sampled n x n stochastic matrix.

##### Arguments

- `n::Integer` : Number of states.
- `k::Integer` : Number of nonzero entries in each row of the matrix.

##### Returns

- `p::Array` : Stochastic matrix.

"""
function random_stochastic_matrix(n::Integer)
    if n <= 0
        throw(ArgumentError("n must be a positive integer"))
    end

    p = random_probvec(n, n)

    return transpose(p)
end

"""
Return a randomly sampled n x n stochastic matrix with k nonzero entries for
each row.

##### Arguments

- `n::Integer` : Number of states.
- `k::Integer` : Number of nonzero entries in each row of the matrix.

##### Returns

- `p::Array` : Stochastic matrix.

"""
function random_stochastic_matrix(n::Integer, k::Integer)
    if !(n > 0)
        throw(ArgumentError("n must be a positive integer"))
    end
    if !(k > 0 && k <= n)
        throw(ArgumentError("k must be an integer with 0 < k <= n"))
    end

    k == n && return random_stochastic_matrix(n)

    # if k < n
    probvecs = random_probvec(k, n)

    # Randomly sample row indices for each column for nonzero values
    row_indices = Array(Int, k*n)
    for j in 1:n
        row_indices[(j-1)*k+1:j*k] = sample(1:n, k, replace=false)
    end

    p = zeros(n, n)
    for j in 1:n
        for i in 1:k
            p[row_indices[(j-1)*k+i], j] = probvecs[i, j]
        end
    end

    return transpose(p)
end

# random_probvec

"""
Return m randomly sampled probability vectors of size k.

##### Arguments

- `k::Integer` : Number of probability vectors.
- `m::Integer` : Size of each probability vectors.

##### Returns

- `a::Array` : Array of shape (k, m) containing probability vectors as colums.

"""
function random_probvec(k::Integer, m::Integer)
    x = Array(Float64, (k+1, m))

    r = rand(k-1, m)
    x[1, :], x[2:end-1, :], x[end, :] = 0, sort(r, 1), 1
    return diff(x, 1)
end


text:
# -*- coding: utf-8 -*-
# testing
import Pkg
Pkg.activate(dirname(Base.current_project())) # activate project
Pkg.precompile()
#Pkg.resolve()
# Pkg.instantiate()

using ODINN
using Test
# using Random
# using Distributed
# using Statistics, LinearAlgebra, Random, Polynomials
# using HDF5
# using JLD2
# using OrdinaryDiffEq, DiffEqFlux
# using Zygote: @ignore
# using Flux
# using Tullio, RecursiveArrayTools
# using Infiltrator
# using Plots
# using ProgressMeter, ParallelDataTransfer
# using Dates
# using PyCall
# using Makie, CairoMakie

# @everywhere begin
#     import Pkg
#     Pkg.activate(dirname(Base.current_project()))
# end

# include all functions
# include("../scripts/toy_model/helpers/parameters.jl")
# include("../test/utils_test.jl")
# include("../scripts/toy_model/helpers/iceflow.jl")

@testset "Halfar" begin

# ######################################################################################################
# ######################             Testing Forward SIA Model            ##############################
# ######################################################################################################

# # I CHANGED ICEFLOW TO REMOVE THE VELOCITIES!!!

# I need to bypass Build_PDE_context. That is, create an object that looks like the output
# of this function and pass it to batch_iceflow_PDE

const nx, ny = 100, 100             # Definition of these constants inside the function? const doesn't work --> global
#t₀ = 1
#t₁ = 5
Δx, Δy = 50, 50
#nx, ny = 100, 100
A₀ = 1e-16
h₀ = 500        # [m]
r₀ = 1000       # [m]
ν = (A₀, h₀, r₀)

B = zeros((nx,ny))

# things from glacier_UDE
rng_seed() = MersenneTwister(123) # random seed
solver = Ralston()
#(@isdefined temp_series) || (const temp_series, norm_temp_series = fake_temp_series(t₁))

temps = [0,-2.0,-3.0,-5.0,-10.0,-12.0,-14.0,-15.0,-20.0]
const temp_series = [repeat([tmp], Int(t₁-t₀)+1) for tmp in temps]

# function SIA_solver_test()

    # this should run with different initializations
    # cosntant were suppose to go here

H₀ = [ halfar_solution(t₀, sqrt((Δx * (i - nx/2))^2 + (Δy * (j - ny/2))^2), ν) for i in 1:nx, j in 1:ny]

H_refs = generate_ref_dataset(temp_series, H₀)   # we need to specify the time!!!

# ################ Non-negativity of solutions ##################

function non_negative_solutions(Hs)
    @testset "Non-negativity of solutions" begin
    for i in 1:length(Hs)
        @test 0.0 <= minimum(Hs[i])
    end
end
end

# ############### Conservation of mass ##################

function mass_conservation(H₀, Hs)
    @testset "Conservation of ice mass under zero mass balance conditions" begin
        for i in 1:length(Hs)
            @test norm(Hs[i],1) ≈ norm(H₀,1) atol=1e-2
        end
    end
end

# ############### Halfar Solutions ##################

function halfar_test(H₀, Hs)
    @testset "Halfar Solution Test" begin
    for i in 1:length(H_refs)
        T_mean = mean(temp_series[i])
        println("Current temperature:", T_mean)
        A₁ = A_fake(T_mean)
        τ₁ = t₀ + (A₀ / A₁) * ( t₁ - t₀ )
        H₁ = [ halfar_solution(τ₁, sqrt((Δx * (i - nx/2))^2 + (Δy * (j - ny/2))^2), ν) for i in 1:nx, j in 1:ny];

        println("Norms")
        println( norm(H₀, 1))
        println( norm(H_refs[i],1))
        println( norm(H₁, 1))

        #display(plot(heatmap(H₀), heatmap(H₁), heatmap(H_refs[i]), heatmap(abs.(H₁ - H_refs[i])) ,layout=(2,2)))
        #display(heatmap!(H_refs[i]))

        abs_diff = maximum(abs.(H₁ .- H_refs[i]))

        println(abs_diff)
        @test abs_diff < 0.1
    end
    end
end

### All the test together
non_negative_solutions(H_refs)
mass_conservation(H₀, H_refs)
halfar_test(H₀, H_refs)

# # Test UA where this is just the identity function

UA = FastChain(
    FastDense(1, 1, x->A_fake.(x))
    )
θ = initial_params(UA)

temp = temps[1]

UA(0.0, θ)

predict_A̅(UA, θ, temp_series[1])

end


text:
const t=2.0

function initial_condition_generator(rng,n,model)
    i = rand(rng,1:2,n)
    X = Vector{Float64}(undef,n)
    for k in 1:n
        x = Inf
        while x>=model.b
            x = -log(rand(rng))
        end
        X[k]=x
    end
    return (X=X,φ=i)
end

const tau = fixed_time(t)

const x = range(0,10,length=10_001)


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct FixedAdmittance <: ElectricLoad
    name::String
    available::Bool
    bus::Bus
    Y::Complex{Float64}
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function FixedAdmittance(name, available, bus, Y, _forecasts=InfrastructureSystems.Forecasts(), )
    FixedAdmittance(name, available, bus, Y, _forecasts, InfrastructureSystemsInternal())
end

function FixedAdmittance(; name, available, bus, Y, _forecasts=InfrastructureSystems.Forecasts(), )
    FixedAdmittance(name, available, bus, Y, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function FixedAdmittance(::Nothing)
    FixedAdmittance(;
        name="init",
        available=false,
        bus=Bus(nothing),
        Y=0.0,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get FixedAdmittance name."""
get_name(value::FixedAdmittance) = value.name
"""Get FixedAdmittance available."""
get_available(value::FixedAdmittance) = value.available
"""Get FixedAdmittance bus."""
get_bus(value::FixedAdmittance) = value.bus
"""Get FixedAdmittance Y."""
get_Y(value::FixedAdmittance) = value.Y
"""Get FixedAdmittance _forecasts."""
get__forecasts(value::FixedAdmittance) = value._forecasts
"""Get FixedAdmittance internal."""
get_internal(value::FixedAdmittance) = value.internal


text:
type BatchResultErrorEntry end
type ChangeMessageVisibilityBatchRequestEntry end
type ChangeMessageVisibilityBatchResultEntry end
type DeleteMessageBatchRequestEntry end
type DeleteMessageBatchResultEntry end
type Message end
type MessageAttributeValue end
type SendMessageBatchRequestEntry end
type SendMessageBatchResultEntry end


text:
# Download Prana artifact
using ArtifactUtils, Pkg.Artifacts
add_artifact!(
    "../Artifacts.toml",
    "Prana",
    "https://github.com/stillyslalom/prana/archive/v0.1.0.tar.gz",
    # lazy=true, # need to use LazyArtifacts for this
    force=true,
    clear=false,
)


text:
module TestOperators
    using NullableArrays
    using Base.Test

    A = rand(10) .+ fill(1, 10)
    B = rand(10) .+ fill(1, 10)
    M = rand(Bool, 10)
    i = rand(1:10)
    j = rand(1:10)
    while j == i
        j = rand(1:10)
    end
    M[i] = true
    M[j] = false
    U = NullableArray(A)
    V = NullableArray(A, M)
    X = NullableArray(B)
    Y = NullableArray(B, M)

    srand(1)
    ensure_neg(x::Unsigned) = -convert(Signed, x)
    ensure_neg(x::Any) = -abs(x)

    # check for fast path (null-safe combinations of operators and types)
    if VERSION >= v"0.5.0-dev" # Some functors are missing on 0.4.0, don't check fast path
        for S in NullableArrays.SafeTypes.types,
            T in NullableArrays.SafeTypes.types
            # mixing signed and unsigned types is unsafe (slow path tested below)
            ((S <: Signed) $ (T <: Signed)) && continue

            u0 = zero(S)
            u1 = one(S)
            u2 = rand(S)

            v0 = zero(T)
            v1 = one(T)
            v2 = rand(T)

            # safe unary operators
            for op in (+, -, ~, abs, abs2, cbrt)
                S <: AbstractFloat && op == (~) && continue

                @test op(Nullable(u0)) === Nullable(op(u0))
                @test op(Nullable(u1)) === Nullable(op(u1))
                @test op(Nullable(u2)) === Nullable(op(u2))
                @test op(Nullable(u0, true)) === Nullable(op(u0), true)
            end

            for u in (u0, u1, u2), v in (v0, v1, v2)
                # safe binary operators: === checks that the fast-path was taken (no branch)
                for op in (+, -, *, /, &, |, >>, <<, >>>,
                           <, >, <=, >=,
                           Base.scalarmin, Base.scalarmax)
                    (T <: AbstractFloat || S <: AbstractFloat) && op in (&, |, >>, <<, >>>) && continue

                    @test op(Nullable(u), Nullable(v)) === Nullable(op(u, v))
                    @test op(Nullable(u, true), Nullable(v, true)) === Nullable(op(u, v), true)
                    @test op(Nullable(u), Nullable(v, true)) === Nullable(op(u, v), true)
                    @test op(Nullable(u, true), Nullable(v)) === Nullable(op(u, v), true)
                end
            end
        end

        @test !Nullable(true) === Nullable(false)
        @test !Nullable(false) === Nullable(true)
        @test !(Nullable(true, true)) === Nullable(false, true)
        @test !(Nullable(false, true)) === Nullable(true, true)
    end

    # test all types and operators (including null-unsafe ones)
    for S in Union{NullableArrays.SafeTypes, BigInt, BigFloat}.types,
        T in Union{NullableArrays.SafeTypes, BigInt, BigFloat}.types
        u0 = zero(S)
        u1 = one(S)
        u2 = S <: Union{BigInt, BigFloat} ? S(rand(Int128)) : rand(S)

        v0 = zero(T)
        v1 = one(T)
        v2 = T <: Union{BigInt, BigFloat} ? T(rand(Int128)) : rand(T)

        (v2 > 5 || v2 < -5) && (v2 = T(5)) # Work around JuliaLang/julia#16989

        # safe unary operators
        for op in (+, -, ~, abs, abs2, cbrt)
            T <: AbstractFloat && op == (~) && continue

            R = Base.promote_op(op, T)
            x = op(Nullable(v0))
            @test isa(x, Nullable{R}) && isequal(x, Nullable(op(v0)))
            x = op(Nullable(v1))
            @test isa(x, Nullable{R}) && isequal(x, Nullable(op(v1)))
            x = op(Nullable(v2))
            @test isa(x, Nullable{R}) && isequal(x, Nullable(op(v2)))
            x = op(Nullable(v0, true))
            @test isa(x, Nullable{R}) && isnull(x)
            x = op(Nullable(v1, true))
            @test isa(x, Nullable{R}) && isnull(x)
            x = op(Nullable(v2, true))
            @test isa(x, Nullable{R}) && isnull(x)
            x = op(Nullable{R}())
            @test isa(x, Nullable{R}) && isnull(x)

            x = op(Nullable())
            @test isa(x, Nullable{Union{}}) && isnull(x)
        end

        # unsafe unary operators
        # sqrt
        @test_throws DomainError sqrt(Nullable(ensure_neg(v1)))
        R = Base.promote_op(sqrt, T)
        x = sqrt(Nullable(v0))
        @test isa(x, Nullable{R}) && isequal(x, Nullable(sqrt(v0)))
        x = sqrt(Nullable(v1))
        @test isa(x, Nullable{R}) && isequal(x, Nullable(sqrt(v1)))
        x = sqrt(Nullable(abs(v2)))
        @test isa(x, Nullable{R}) && isequal(x, Nullable(sqrt(abs(v2))))
        x = sqrt(Nullable(v0, true))
        @test isa(x, Nullable{R}) && isnull(x)
        x = sqrt(Nullable(ensure_neg(v1), true))
        @test isa(x, Nullable{R}) && isnull(x)
        x = sqrt(Nullable(ensure_neg(v2), true))
        @test isa(x, Nullable{R}) && isnull(x)
        x = sqrt(Nullable{R}())
        @test isa(x, Nullable{R}) && isnull(x)

        x = sqrt(Nullable())
        @test isa(x, Nullable{Union{}}) && isnull(x)

        for u in (u0, u1, u2), v in (v0, v1, v2)
            # safe binary operators
            for op in (+, -, *, /, &, |, >>, <<, >>>,
                       <, >, <=, >=,
                       Base.scalarmin, Base.scalarmax)
                (T <: AbstractFloat || S <: AbstractFloat) && op in (&, |, >>, <<, >>>) && continue
                if VERSION < v"0.5.0-dev"
                    (T <: Bool || S <: Bool) && op in (>>, <<, >>>) && continue
                    (T <: BigInt || S <: BigInt) && op in (&, |, >>, <<, >>>) && continue
                end

                if S <: Unsigned || T <: Unsigned
                    @test isequal(op(Nullable(abs(u)), Nullable(abs(v))), Nullable(op(abs(u), abs(v))))
                else
                    @test isequal(op(Nullable(u), Nullable(v)), Nullable(op(u, v)))
                end
                R = Base.promote_op(op, S, T)
                x = op(Nullable(u, true), Nullable(v, true))
                @test isa(x, Nullable{R}) && isnull(x)
                x = op(Nullable(u), Nullable(v, true))
                @test isa(x, Nullable{R}) && isnull(x)
                x = op(Nullable(u, true), Nullable(v))
                @test isa(x, Nullable{R}) && isnull(x)

                x = op(Nullable(u, true), Nullable())
                @test isa(x, Nullable{S}) && isnull(x)
                x = op(Nullable(), Nullable(u, true))
                @test isa(x, Nullable{S}) && isnull(x)
                x = op(Nullable(), Nullable())
                @test isa(x, Nullable{Union{}}) && isnull(x)
            end

            # unsafe binary operators
            # ^
            if S <: Integer && T <: Integer && u != 0 && u != 1 && v != 0
                @test_throws DomainError Nullable(u)^Nullable(ensure_neg(v))
            end
            @test isequal(Nullable(u)^Nullable(one(T)+one(T)), Nullable(u^(one(T)+one(T))))
            R = Base.promote_op(^, S, T)
            x = Nullable(u, true)^Nullable(-abs(v), true)
            @test isnull(x) && eltype(x) === R
            x = Nullable(u, false)^Nullable(-abs(v), true)
            @test isnull(x) && eltype(x) === R
            x = Nullable(u, true)^Nullable(-abs(v), false)
            @test isnull(x) && eltype(x) === R

            x = Nullable(u, true)^Nullable()
            @test isa(x, Nullable{S}) && isnull(x)
            x = Nullable()^Nullable(u, true)
            @test isa(x, Nullable{S}) && isnull(x)
            x = Nullable()^Nullable()
            @test isa(x, Nullable{Union{}}) && isnull(x)

            # ÷ and %
            for op in (÷, %)
                if S <: Integer && T <: Integer && v == 0
                    @test_throws DivideError op(Nullable(u), Nullable(v))
                else
                    @test isequal(op(Nullable(u), Nullable(v)), Nullable(op(u, v)))
                end
                R = Base.promote_op(op, S, T)
                x = op(Nullable(u, true), Nullable(v, true))
                @test isnull(x) && eltype(x) === R
                x = op(Nullable(u, false), Nullable(v, true))
                @test isnull(x) && eltype(x) === R
                x = op(Nullable(u, true), Nullable(v, false))
                @test isnull(x) && eltype(x) === R

                x = op(Nullable(u, true), Nullable())
                @test isa(x, Nullable{S}) && isnull(x)
                x = op(Nullable(), Nullable(u, true))
                @test isa(x, Nullable{S}) && isnull(x)
                x = op(Nullable(), Nullable())
                @test isa(x, Nullable{Union{}}) && isnull(x)
            end

            # isless
            @test isless(Nullable(u), Nullable(v)) === isless(u, v)

            @test isless(Nullable(u), Nullable(v, true)) === true
            @test isless(Nullable(u, true), Nullable(v)) === false
            @test isless(Nullable(u, true), Nullable(v, true)) === false

            @test isless(Nullable(u), Nullable()) === true
            @test isless(Nullable(), Nullable(v)) === false
            @test isless(Nullable(), Nullable()) === false
        end
    end
end


text:
# include("../instance.jl")
# include("../solution.jl")

# module LpTiming
# Déclaration des packages utilisés dans ce fichier
# certains sont déjà chargés dans le fichier usings.jl

# import JuMP
# using .Instance: Instance
# using .Solution: Solution, solve_to_earliest
# using .Model: new_model

# export LpTimingSolver

"""
Ce solveur résoud le  sous-problème de timing consistant à trouver les dates
optimales d'atterrissage des avions à ordre fixé.
Par rapport aux autres solvers (e.g DescentSolver, AnnealingSolver, ...), il
ne contient pas d'attribut bestsol
"""
mutable struct LpTimingSolver
    inst::Instance
    loglevel::Int         # niveau de verbosité
    # Les attributs spécifiques au modèle
    model::Model  # Le modèle MIP
    x         # vecteur des variables d'atterrissage
    cost      # variable du coût de la solution
    costs     # variables du coût de chaque avion

    nb_calls::Int  # POUR FAIRE VOS MESURES DE PERFORMANCE !
    nb_infeasable::Int

    # Le constructeur
    function LpTimingSolver(inst::Instance)
        model = new_model()
        return new(inst, 0, model, zeros(Int, inst.nb_planes),
                   0., zeros(Float64, inst.nb_planes), 0, 0)
    end
end

"""Permet de retrouver le nom de notre XxxxTimingSolver à partir de l'objet"""
function symbol(sv::LpTimingSolver)
    return :lp
end

"""
Trouve les temps d'atterrissage optimales pour un ordre d'avions donné par
la solution `sol`.

Les temps d'atterrissage de la solution sont modifiés in-place.

Cette fonction résoud un problème linéaire. Si le problème n'a pas de solution,
la stratégie d'atterrissage au plus tôt est appliquée.
"""
function solve!(sv::LpTimingSolver, sol::Solution)
    sv.nb_calls += 1
    model = new_model()

    n = sv.inst.nb_planes

    # 1. Création du modèle spécifiquement pour cet ordre d'avion de cette solution
    @variable(model, x[1:n] >= 0)
    @variable(model, y[1:n] >= 0)
    @variable(model, z[1:n] >= 0)

    @objective(model, Min,
               sum(plane.ep * y[i] + plane.tp * z[i]
                   for (i, plane) in enumerate(sol.planes)))

    for (i, plane) in enumerate(sol.planes)
        # Landing time bounds
        @constraint(model, x[i] >= plane.lb)
        @constraint(model, x[i] <= plane.ub)
        # Linearization
        @constraint(model, y[i] >= plane.target - x[i])
        @constraint(model, z[i] >= x[i] - plane.target)
        # Separation with planes landing later
        for (j, second_plane) in enumerate(sol.planes[i+1:end])
            @constraint(model, x[i+j] >= x[i] + get_sep(sv.inst, plane,
                                                        second_plane))
        end
    end

    # 2. résolution du problème à permu d'avion fixée
    JuMP.optimize!(model)

    # 3. Test de la validité du résultat
    if  JuMP.termination_status(model) == MOI.OPTIMAL
        # tout va bien, on peut exploiter le résultat

        # 4. Extraction des valeurs des variables d'atterrissage
        for (i, plane) in enumerate(sol.planes)
            sol.x[i] = round(Int, value(x[i]))
            sol.costs[i] = plane.ep * value(y[i]) + plane.tp * value(z[i])
        end
        prec = Args.args[:cost_precision]
        sol.cost = round(objective_value(model), digits=prec)
    else
        # La solution du solver est invalide : on utilise le placement au plus
        # tôt de façon à disposer malgré tout d'un coût pénalisé afin de pouvoir
        # continuer la recherche heuristique de solutions.
        sv.nb_infeasable += 1
        solve_to_earliest!(sol)
    end
end

# end


text:
function setup_etc(backend::GPUBackend, layer::AccuracyLayer, op_dim::Int, inputs)
  dims = [size(inputs[1])...]
  dims[op_dim] = 1
  etc = make_blob(backend, eltype(inputs[1]), dims...)
  return etc
end
function shutdown(backend::GPUBackend, state::AccuracyLayerState)
  destroy(state.etc)
end

function forward(backend::GPUBackend, state::AccuracyLayerState, inputs::Vector{Blob})
  pred = inputs[1]
  label = inputs[2]

  spatial_dim, pred_dim, num = split_dims(pred, state.op_dim)
  data_type = eltype(pred)

  x_block = round(Int, ceil(convert(Float64, num)/CUDA.THREADS_PER_BLOCK_X));
  y_block = round(Int, ceil(convert(Float64, spatial_dim)/CUDA.THREADS_PER_BLOCK_Y));

  if data_type == Float32
    kernel = backend.mocha.accuracy_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.accuracy_forward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, (x_block,y_block),(CUDA.THREADS_PER_BLOCK_X,CUDA.THREADS_PER_BLOCK_Y),
      (pred.ptr.p, label.ptr.p, state.etc.ptr.p, num, pred_dim, spatial_dim));

  N = num * spatial_dim
  accuracy = CuBLAS.dot(backend.cublas_ctx, data_type, N, state.etc.ptr, 1, state.etc.ptr, 1)

  # accumulate accuracy
  state.accuracy = (state.accuracy * state.n_accum + accuracy) / (N + state.n_accum)
  state.n_accum += N
end


text:
import Dance

include("./utils/main.jl")

function _build_path(sub_dir::String="") :: String
    path::String = abspath(pwd())

    if length(sub_dir) > 0
        path = joinpath(abspath(pwd()), sub_dir)
    end

    if Sys.iswindows()
        path = replace(path, "/" => "\\")
    end

    return path
end

# Test project settings check
@testset "Settings - Project Settings Check" begin
    # Default Prod env does not have :server_host & :server_port defined
    Dance.start_project("demo")
    @test_logs (:error, "Please set valid server_host value, cannot be blank") (:error, "Please set valid server_port value, cannot be blank") Dance.populate_settings(joinpath(abspath(@__DIR__), "demo"))

    project_settings()
    Dance.populate_settings(joinpath(abspath(@__DIR__), "demo"))
    @test haskey(Dance.Configuration.Settings, :foo)
    @test Dance.Configuration.Settings[:foo]=="bar"

    delete_project(false)
end

# Test LOAD_PATH: single-dir project case
@testset "Settings - LOAD_PATH single-dir" begin
    Dance.start_project("demo")
    project_settings()
    dirs_add_single()

    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(5)

    @test _build_path() in LOAD_PATH
    @test _build_path("html") in LOAD_PATH
    @test _build_path("settings") in LOAD_PATH
    @test _build_path("code") in LOAD_PATH

    delete_project()
end

# Test LOAD_PATH: multiple sub-dirs case
@testset "Settings - LOAD_PATH multiple sub-dirs" begin
    Dance.start_project("demo")
    project_settings()
    dirs_add_multiple()

    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(2)

    @test _build_path() in LOAD_PATH
    @test _build_path("html") in LOAD_PATH
    @test _build_path("settings") in LOAD_PATH
    @test _build_path("code") in LOAD_PATH
    @test _build_path("code/sub-dir1") in LOAD_PATH
    @test _build_path("code/sub-dir1/sub-dir2") in LOAD_PATH

    delete_project()
end

# Test `files/routes.jl` default contents
@testset "Settings - default routes contents" begin
    Dance.start_project("demo")
    project_settings()

    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(1)

    make_and_test_request_get("/", 200, Dict("Content-Type" => "text/html; charset=UTF-8"), 217, false, "Hello World")

    delete_project()
end


text:
module Filter

using DSP, FFTW
using Match

export decode, filter_signal, get_filter

get_symbol(f1, f2) = @match (f1, f2) begin

	(941, 1336) => "0"
	(697, 1209) => "1"
	(697, 1336) => "2"
	(697, 1477) => "3"
	(770, 1209) => "4"
	(770, 1336) => "5"
	(770, 1477) => "6"
	(852, 1209) => "7"
	(852, 1336) => "8"
	(852, 1477) => "9"
	(941, 1209) => "*"
	(941, 1477) => "#"
	(697, 1633) => "A"
	(770, 1633) => "B"
	(852, 1633) => "C"
	(941, 1633) => "D"
	_ => ""
end

function get_filter(fs, width, fc)

	x_h = width .* (-0.05:(1/fs):0.05)

	h = sinc.(x_h)

	wc = ((2*pi)/fs) * fc

	h_bp = h .* cos.((0:(length(h) - 1)) .* wc)

	return h_bp
end

# Filter the signal for
# the given frecuency 'fc'
function filter_signal(s, fs; width, fc)

	# Create the filter
	bp = get_filter(fs, width, fc)

	# Filter the signal
	y = conv(s, bp)

	return y
end

function get_window(fs, f)

	N = Int(round((1/f) * fs)) * 3;

	win = tukey(N, 0.50)

	return win
end

function get_envolvent(win, e)

	# Envolvent
	y = conv(win, e)

	y = y[div(length(win), 2):end]

	# Normalize envolvent
	y = y .* (maximum(abs.(e))/maximum(abs.(y)))

	return y
end

function detect_pulses(y)

	# Threshold
	thres = ones(length(y)) .* (1.45)

	# Detect frequency pulses from the other things
	pulses = map(x -> Int(x), y .> thres)

	# Detect where the pulses
	# start and end
	differ = diff(pulses)

	# Indexes where the pulses start
	pulses_start = map(x -> x + 1, findall(x -> x == 1, differ))

	return pulses_start
end

function analyze_pulses(lowy, highy, lowf, highf, fs)

	# Energies
	lowe = lowy .^ 2
	highe = highy .^ 2

	# Store max in energy to
	# compare at the end
	maxe = (maximum(abs.(lowe)) + maximum(abs.(highe))) / 2

	# Normalize energies
	lowe = lowe ./ maximum(abs.(lowe))
	highe = highe ./ maximum(abs.(highe))

	# Sum of the energies
	e = lowe .+ highe

	win = get_window(fs, lowf)

	y = get_envolvent(win, e)

	pulses = detect_pulses(y)

	# Symbol mapping
	sym = get_symbol(lowf, highf)
	pulses = map(x -> (sym, x, maxe), pulses)

	return pulses
end

function analyze_energy(symbols)

	# Analize maximums of energy
	max_energies = map(x -> x[3], symbols)
	maxe = maximum(max_energies)

	# Normalize energy maximums
	symbols = map(x -> (x[1], x[2], x[3]/maxe), symbols)

	# Discard the elements that do not pass
	# the energy threshold of 20%
	symbols = filter(x -> x[3] > 0.20, symbols)

	# Sort by sample starting
	sort!(symbols, by = x -> x[2])

	symbols = map(x -> x[1], symbols)

	return symbols
end

function decode(s, fs)

	s = vec(s)

	# Combine low and high frequencies
	# and analiyze them to decode the symbols
	#
	# After this analysis, 'symbols' will
	# contain tuples in the following form:
	#
	# syms = [(3, 500), (2,1000), (3,400), (3,800), (7,1200)]
	#
	# and then it will be sorted by the second element
	# of the tuple, which is the sample where it begins that
	# pulse
	#
	# syms = [(3,400), (2,500), (3,800), (2,1000), (7,1200)]
	#
	# Finally a 'map' operation will remove the
	# samples startings,
	#
	# syms = [3, 2, 3, 2, 7]
	#

	symbols = []

	high_frecs = Dict([f => filter_signal(s, fs, fc=f, width=50)
			   for f in [1209, 1336, 1477, 1633]])

	for lowf in [697, 770, 852, 941]

		lowy = filter_signal(s, fs, fc=lowf, width=40)

		for highf in [1209, 1336, 1477, 1633]

			highy = high_frecs[highf]

			syms = analyze_pulses(lowy, highy, lowf, highf, fs)

			append!(symbols, syms)
		end
	end

	symbols = analyze_energy(symbols)

	return symbols
end

end # module

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
"""
    size(A) -> Tuple
    size(A, dim) -> Union{Int,StaticInt}

Returns the size of each dimension of `A` or along dimension `dim` of `A`. If the size of
any axes are known at compile time, these should be returned as `Static` numbers. Otherwise,
`ArrayInterface.size(A)` is identical to `Base.size(A)`

```julia
julia> using StaticArrays, ArrayInterface

julia> A = @SMatrix rand(3,4);

julia> ArrayInterface.size(A)
(static(3), static(4))
```
"""
size(a::A) where {A} = _maybe_size(Base.IteratorSize(A), a)
size(a::Base.Broadcast.Broadcasted) = map(length, axes(a))

_maybe_size(::Base.HasShape{N}, a::A) where {N,A} = map(length, axes(a))
_maybe_size(::Base.HasLength, a::A) where {A} = (length(a),)
size(x::SubArray) = eachop(_sub_size, to_parent_dims(x), x.indices)
_sub_size(x::Tuple, ::StaticInt{dim}) where {dim} = length(getfield(x, dim))
@inline size(B::VecAdjTrans) = (One(), length(parent(B)))
@inline size(B::MatAdjTrans) = permute(size(parent(B)), to_parent_dims(B))
@inline function size(B::PermutedDimsArray{T,N,I1}) where {T,N,I1}
    permute(size(parent(B)), static(I1))
end
function size(a::ReinterpretArray{T,N,S,A}) where {T,N,S,A}
    psize = size(parent(a))
    if _is_reshaped(typeof(a))
        if sizeof(S) === sizeof(T)
            return psize
        elseif sizeof(S) > sizeof(T)
            return (static(div(sizeof(S), sizeof(T))), psize...)
        else
            return tail(psize)
        end
    else
        return (div(first(psize) * static(sizeof(S)), static(sizeof(T))), tail(psize)...,)
    end
end
size(A::ReshapedArray) = Base.size(A)
size(A::AbstractRange) = (length(A),)
size(x::Base.Generator) = size(getfield(x, :iter))
size(x::Iterators.Reverse) = size(getfield(x, :itr))
size(x::Iterators.Enumerate) = size(getfield(x, :itr))
size(x::Iterators.Accumulate) = size(getfield(x, :itr))
size(x::Iterators.Pairs) = size(getfield(x, :itr))
@inline function size(x::Iterators.ProductIterator)
    eachop(_sub_size, nstatic(Val(ndims(x))), getfield(x, :iterators))
end

size(a, dim) = size(a, to_dims(a, dim))
size(a::Array, dim::Integer) = Base.arraysize(a, convert(Int, dim))
function size(a::A, dim::Integer) where {A}
    if parent_type(A) <: A
        len = known_size(A, dim)
        if len === nothing
            return Int(length(axes(a, dim)))
        else
            return StaticInt(len)
        end
    else
        return size(a)[dim]
    end
end
function size(A::SubArray, dim::Integer)
    pdim = to_parent_dims(A, dim)
    if pdim > ndims(parent_type(A))
        return size(parent(A), pdim)
    else
        return length(A.indices[pdim])
    end
end
size(x::Iterators.Zip) = Static.reduce_tup(promote_shape, map(size, getfield(x, :is)))

"""
    known_size(::Type{T}) -> Tuple
    known_size(::Type{T}, dim) -> Union{Int,Nothing}

Returns the size of each dimension of `A` or along dimension `dim` of `A` that is known at
compile time. If a dimension does not have a known size along a dimension then `nothing` is
returned in its position.
"""
known_size(x) = known_size(typeof(x))
function known_size(::Type{T}) where {T<:AbstractRange}
    (_range_length(known_first(T), known_step(T), known_last(T)),)
end
known_size(::Type{<:Base.Generator{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Reverse{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Enumerate{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Accumulate{<:Any,I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Pairs{<:Any,<:Any,I}}) where {I} = known_size(I)
@inline function known_size(::Type{<:Iterators.ProductIterator{T}}) where {T}
    eachop(_known_size, nstatic(Val(known_length(T))), T)
end

# 1. `Zip` doesn't check that its collections are compatible (same size) at construction,
#   but we assume as much b/c otherwise it will error while iterating. So we promote to the
#   known size if matching a `Nothing` and `Int` size.
# 2. `promote_shape(::Tuple{Vararg{CanonicalInt}}, ::Tuple{Vararg{CanonicalInt}})` promotes
#   trailing dimensions (which must be of size 1), to `static(1)`. We want to stick to
#   `Nothing` and `Int` types, so we do one last pass to ensure everything is dynamic
@inline function known_size(::Type{<:Iterators.Zip{T}}) where {T}
    dynamic(reduce_tup(_promote_shape, eachop(_unzip_size, nstatic(Val(known_length(T))), T)))
end
_unzip_size(::Type{T}, n::StaticInt{N}) where {T,N} = known_size(field_type(T, n))

known_size(::Type{T}) where {T} = _maybe_known_size(Base.IteratorSize(T), T)
function _maybe_known_size(::Base.HasShape{N}, ::Type{T}) where {N,T}
    eachop(_known_size, nstatic(Val(N)), axes_types(T))
end
_maybe_known_size(::Base.IteratorSize, ::Type{T}) where {T} = (known_length(T),)
_known_size(::Type{T}, dim::StaticInt) where {T} = known_length(field_type(T, dim))
@inline known_size(x, dim) = known_size(typeof(x), dim)
@inline known_size(::Type{T}, dim) where {T} = known_size(T, to_dims(T, dim))
@inline function known_size(::Type{T}, dim::CanonicalInt) where {T}
    if ndims(T) < dim
        return 1
    else
        return known_size(T)[dim]
    end
end


text:
"""
    save_fit(bbl::Union{GlobalModel, BlackBoxLearner, Array}, dir::String)

Saves IAI fits associated with different OptimalConstraintTree objects.
"""
function save_fit(bbl::BlackBoxLearner, dir::String = TREE_DIR)
    # TODO: support saving all of the BBR learner history.
    IAI.write_json(dir * bbl.name * ".json", bbl.learners[end]) # save learner
    if bbl isa BlackBoxClassifier
        save(dir * bbl.name * ".jld", Dict("learner_kwargs" => bbl.learner_kwargs[end],
                                        "accuracies" => bbl.accuracies[end]))
    else
        save(dir * bbl.name * ".jld", Dict("learner_kwargs" => bbl.learner_kwargs[end],
                                        "thresholds" => bbl.thresholds[end],
                                        "ul_data" => bbl.ul_data[end]))
    end
    return
end

save_fit(bbls::Array, dir::String = TREE_DIR) = [save_fit(bbl, dir) for bbl in bbls]
save_fit(gm::GlobalModel, dir::String = TREE_DIR) = save_fit(gm.bbls, dir * "$(gm.name)_")

"""
    load_fit(BlackBoxLearner, dir::String = TREE_DIR)

Loads IAI fits associated with OptimalConstraintTree objects.
Checks that there is correspondence between loaded trees and the associated constraints.
"""

function load_fit(bbl::BlackBoxLearner, dir::String = TREE_DIR)
    loaded_lnr = IAI.read_json(dir * bbl.name * ".json");
    size(IAI.variable_importance(loaded_lnr), 1) == length(bbl.vars) || throw(
        OCTException("Object " * bbl.name * " does not match associated learner."))
    set_param(bbl, :reloaded, true)
    push!(bbl.learners, loaded_lnr)
    dd = load(dir * bbl.name * ".jld")
    if bbl isa BlackBoxClassifier
        push!(bbl.learner_kwargs, dd["learner_kwargs"])
        push!(bbl.accuracies, dd["accuracies"])
    else
        push!(bbl.learner_kwargs, dd["learner_kwargs"])
        push!(bbl.thresholds, dd["thresholds"])
        push!(bbl.ul_data, dd["ul_data"])
    end
    return
end

load_fit(bbls::Array{BlackBoxLearner},
        dir::String = TREE_DIR) = [load_fit(bbl, dir) for bbl in bbls]

load_fit(gm::GlobalModel, dir::String = TREE_DIR) = load_fit(gm.bbls, dir * "$(gm.name)_")


text:
@hascuda using CUDAnative, CuArrays

import GPUifyLoops: @launch, @loop, @unroll, @synchronize

using Oceananigans.Operators

const Tx = 16 # CUDA threads per x-block
const Ty = 16 # CUDA threads per y-block

"""
    time_step!(model, Nt, Δt)

Step forward `model` `Nt` time steps using a second-order Adams-Bashforth
method with step size `Δt`.
"""
function time_step!(model::Model{A}, Nt, Δt) where A <: Architecture
    clock = model.clock
    model_start_time = clock.time
    model_end_time = model_start_time + Nt*Δt

    if clock.iteration == 0
        for output_writer in model.output_writers
            write_output(model, output_writer)
        end
        for diagnostic in model.diagnostics
            run_diagnostic(model, diagnostic)
        end
    end

    arch = A()

    Nx, Ny, Nz = model.grid.Nx, model.grid.Ny, model.grid.Nz
    Lx, Ly, Lz = model.grid.Lx, model.grid.Ly, model.grid.Lz
    Δx, Δy, Δz = model.grid.Δx, model.grid.Δy, model.grid.Δz

    grid = model.grid
    cfg = model.configuration
    bcs = model.boundary_conditions
    clock = model.clock

    G = model.G
    Gp = model.Gp
    constants = model.constants
    eos =  model.eos
    U = model.velocities
    tr = model.tracers
    pr = model.pressures
    forcing = model.forcing
    poisson_solver = model.poisson_solver

    RHS = model.stepper_tmp.fCC1
    ϕ = model.stepper_tmp.fCC2

    gΔz = model.constants.g * model.grid.Δz
    fCor = model.constants.f

    uvw = U.u.data, U.v.data, U.w.data
    TS = tr.T.data, tr.S.data
    Guvw = G.Gu.data, G.Gv.data, G.Gw.data

    # Source terms at current (Gⁿ) and previous (G⁻) time steps.
    Gⁿ = G.Gu.data, G.Gv.data, G.Gw.data, G.GT.data, G.GS.data
    G⁻ = Gp.Gu.data, Gp.Gv.data, Gp.Gw.data, Gp.GT.data, Gp.GS.data

    Bx, By, Bz = floor(Int, Nx/Tx), floor(Int, Ny/Ty), Nz  # Blocks in grid

    tb = (threads=(Tx, Ty), blocks=(Bx, By, Bz))

    for n in 1:Nt
        χ = ifelse(model.clock.iteration == 0, -0.5, 0.125) # Adams-Bashforth (AB2) parameter.

        @launch device(arch) store_previous_source_terms!(grid, Gⁿ..., G⁻..., threads=(Tx, Ty), blocks=(Bx, By, Bz))
        @launch device(arch) update_buoyancy!(grid, constants, eos, tr.T.data, pr.pHY′.data, threads=(Tx, Ty), blocks=(Bx, By))
        @launch device(arch) calculate_interior_source_terms!(grid, constants, eos, cfg, uvw..., TS..., pr.pHY′.data, Gⁿ..., forcing, threads=(Tx, Ty), blocks=(Bx, By, Bz))
                             calculate_boundary_source_terms!(model)
        @launch device(arch) adams_bashforth_update_source_terms!(grid, Gⁿ..., G⁻..., χ, threads=(Tx, Ty), blocks=(Bx, By, Bz))
        @launch device(arch) calculate_source_term_divergence!(arch, grid, Guvw..., RHS.data, threads=(Tx, Ty), blocks=(Bx, By, Bz))

        if arch == CPU()
            solve_poisson_3d_ppn_planned!(poisson_solver, grid, RHS, ϕ)
            @. pr.pNHS.data = real(ϕ.data)
        elseif arch == GPU()
            solve_poisson_3d_ppn_gpu_planned!(Tx, Ty, Bx, By, Bz, poisson_solver, grid, RHS, ϕ)
            @launch device(arch) idct_permute!(grid, ϕ.data, pr.pNHS.data, threads=(Tx, Ty), blocks=(Bx, By, Bz))
        end

        @launch device(arch) update_velocities_and_tracers!(grid, uvw..., TS..., pr.pNHS.data, Gⁿ..., G⁻..., Δt, threads=(Tx, Ty), blocks=(Bx, By, Bz))

        clock.time += Δt
        clock.iteration += 1

        for diagnostic in model.diagnostics
            (clock.iteration % diagnostic.diagnostic_frequency) == 0 && run_diagnostic(model, diagnostic)
        end

        for output_writer in model.output_writers
            (clock.iteration % output_writer.output_frequency) == 0 && write_output(model, output_writer)
        end
    end

    return nothing
end

time_step!(model; Nt, Δt) = time_step!(model, Nt, Δt)

"""Store previous source terms before updating them."""
function store_previous_source_terms!(grid::Grid, Gu, Gv, Gw, GT, GS, Gpu, Gpv, Gpw, GpT, GpS)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gpu[i, j, k] = Gu[i, j, k]
                @inbounds Gpv[i, j, k] = Gv[i, j, k]
                @inbounds Gpw[i, j, k] = Gw[i, j, k]
                @inbounds GpT[i, j, k] = GT[i, j, k]
                @inbounds GpS[i, j, k] = GS[i, j, k]
            end
        end
    end
    @synchronize
end

"Update the hydrostatic pressure perturbation pHY′ and buoyancy δρ."
function update_buoyancy!(grid::Grid, constants, eos, T, pHY′)
    gΔz = constants.g * grid.Δz

    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            @inbounds pHY′[i, j, 1] = 0.5 * gΔz * δρ(eos, T, i, j, 1)
            @unroll for k in 2:grid.Nz
                @inbounds pHY′[i, j, k] = pHY′[i, j, k-1] + gΔz * 0.5 * (δρ(eos, T, i, j, k-1) + δρ(eos, T, i, j, k))
            end
        end
    end

    @synchronize
end

"Store previous value of the source term and calculate current source term."
function calculate_interior_source_terms!(grid::Grid, constants, eos, cfg, u, v, w, T, S, pHY′, Gu, Gv, Gw, GT, GS, F)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    Δx, Δy, Δz = grid.Δx, grid.Δy, grid.Δz

    fCor = constants.f
    ρ₀ = eos.ρ₀
    𝜈h, 𝜈v, κh, κv = cfg.𝜈h, cfg.𝜈v, cfg.κh, cfg.κv

    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                # u-momentum equation
                @inbounds Gu[i, j, k] = (-u∇u(grid, u, v, w, i, j, k)
                                            + fCor*avg_xy(grid, v, i, j, k)
                                            - δx_c2f(grid, pHY′, i, j, k) / (Δx * ρ₀)
                                            + 𝜈∇²u(grid, u, 𝜈h, 𝜈v, i, j, k)
                                            + F.u(grid, u, v, w, T, S, i, j, k))

                # v-momentum equation
                @inbounds Gv[i, j, k] = (-u∇v(grid, u, v, w, i, j, k)
                                            - fCor*avg_xy(grid, u, i, j, k)
                                            - δy_c2f(grid, pHY′, i, j, k) / (Δy * ρ₀)
                                            + 𝜈∇²v(grid, v, 𝜈h, 𝜈v, i, j, k)
                                            + F.v(grid, u, v, w, T, S, i, j, k))

                # w-momentum equation: comment about how pressure and buoyancy are handled
                @inbounds Gw[i, j, k] = (-u∇w(grid, u, v, w, i, j, k)
                                            + 𝜈∇²w(grid, w, 𝜈h, 𝜈v, i, j, k)
                                            + F.w(grid, u, v, w, T, S, i, j, k))

                # temperature equation
                @inbounds GT[i, j, k] = (-div_flux(grid, u, v, w, T, i, j, k)
                                            + κ∇²(grid, T, κh, κv, i, j, k)
                                            + F.T(grid, u, v, w, T, S, i, j, k))

                # salinity equation
                @inbounds GS[i, j, k] = (-div_flux(grid, u, v, w, S, i, j, k)
                                            + κ∇²(grid, S, κh, κv, i, j, k)
                                            + F.S(grid, u, v, w, T, S, i, j, k))
            end
        end
    end

    @synchronize
end

function adams_bashforth_update_source_terms!(grid::Grid, Gu, Gv, Gw, GT, GS, Gpu, Gpv, Gpw, GpT, GpS, χ)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gu[i, j, k] = (1.5 + χ)*Gu[i, j, k] - (0.5 + χ)*Gpu[i, j, k]
                @inbounds Gv[i, j, k] = (1.5 + χ)*Gv[i, j, k] - (0.5 + χ)*Gpv[i, j, k]
                @inbounds Gw[i, j, k] = (1.5 + χ)*Gw[i, j, k] - (0.5 + χ)*Gpw[i, j, k]
                @inbounds GT[i, j, k] = (1.5 + χ)*GT[i, j, k] - (0.5 + χ)*GpT[i, j, k]
                @inbounds GS[i, j, k] = (1.5 + χ)*GS[i, j, k] - (0.5 + χ)*GpS[i, j, k]
            end
        end
    end
    @synchronize
end

"Store previous value of the source term and calculate current source term."
function calculate_source_term_divergence!(::CPU, grid::Grid, Gu, Gv, Gw, RHS)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                # Calculate divergence of the RHS source terms (Gu, Gv, Gw).
                @inbounds RHS[i, j, k] = div_f2c(grid, Gu, Gv, Gw, i, j, k)
            end
        end
    end

    @synchronize
end

function calculate_source_term_divergence!(::GPU, grid::Grid, Gu, Gv, Gw, RHS)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:Nz; blockIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                # Calculate divergence of the RHS source terms (Gu, Gv, Gw) and applying a permutation
                # which is the first step in the DCT.
                if CUDAnative.ffs(k) == 1  # isodd(k)
                    @inbounds RHS[i, j, convert(UInt32, CUDAnative.floor(k/2) + 1)] = div_f2c(grid, Gu, Gv, Gw, i, j, k)
                else
                    @inbounds RHS[i, j, convert(UInt32, Nz - CUDAnative.floor((k-1)/2))] = div_f2c(grid, Gu, Gv, Gw, i, j, k)
                end
            end
        end
    end

    @synchronize
end

function idct_permute!(grid::Grid, ϕ, pNHS)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:Nz; blockIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                if k <= Nz/2
                    @inbounds pNHS[i, j, 2k-1] = real(ϕ[i, j, k])
                else
                    @inbounds pNHS[i, j, 2(Nz-k+1)] = real(ϕ[i, j, k])
                end
            end
        end
    end

    @synchronize
end

function update_velocities_and_tracers!(grid::Grid, u, v, w, T, S, pNHS, Gu, Gv, Gw, GT, GS, Gpu, Gpv, Gpw, GpT, GpS, Δt)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds u[i, j, k] = u[i, j, k] + (Gu[i, j, k] - (δx_c2f(grid, pNHS, i, j, k) / grid.Δx)) * Δt
                @inbounds v[i, j, k] = v[i, j, k] + (Gv[i, j, k] - (δy_c2f(grid, pNHS, i, j, k) / grid.Δy)) * Δt
                @inbounds w[i, j, k] = w[i, j, k] + (Gw[i, j, k] - (δz_c2f(grid, pNHS, i, j, k) / grid.Δz)) * Δt
                @inbounds T[i, j, k] = T[i, j, k] + (GT[i, j, k] * Δt)
                @inbounds S[i, j, k] = S[i, j, k] + (GS[i, j, k] * Δt)
            end
        end
    end

    @synchronize
end

#
# Boundary condition physics specification
#

"Apply boundary conditions by modifying the source term G."
function calculate_boundary_source_terms!(model::Model{A}) where A <: Architecture
    arch = A()

    Nx, Ny, Nz = model.grid.Nx, model.grid.Ny, model.grid.Nz
    Lx, Ly, Lz = model.grid.Lx, model.grid.Ly, model.grid.Lz
    Δx, Δy, Δz = model.grid.Δx, model.grid.Δy, model.grid.Δz

    grid = model.grid
    clock = model.clock
    eos =  model.eos
    cfg = model.configuration
    bcs = model.boundary_conditions
    U = model.velocities
    tr = model.tracers
    G = model.G

    t, iteration = clock.time, clock.iteration
    u, v, w, T, S = U.u.data, U.v.data, U.w.data, tr.T.data, tr.S.data
    Gu, Gv, Gw, GT, GS = G.Gu.data, G.Gv.data, G.Gw.data, G.GT.data, G.GS.data

    Bx, By, Bz = floor(Int, Nx/Tx), floor(Int, Ny/Ty), Nz  # Blocks in grid

    coord = :z #for coord in (:x, :y, :z) when we are ready to support more coordinates.
    𝜈 = cfg.𝜈v
    κ = cfg.κv

    u_x_bcs = getproperty(bcs.u, coord)
    v_x_bcs = getproperty(bcs.v, coord)
    w_x_bcs = getproperty(bcs.w, coord)
    T_x_bcs = getproperty(bcs.T, coord)
    S_x_bcs = getproperty(bcs.S, coord)

    # Apply boundary conditions. We assume there is one molecular 'diffusivity'
    # value, which is passed to apply_bcs.
    apply_bcs!(arch, Val(coord), Bx, By, Bz, u_x_bcs.left, u_x_bcs.right, u, Gu, 𝜈, u, v, w, T, S, t, iteration, grid)
    apply_bcs!(arch, Val(coord), Bx, By, Bz, v_x_bcs.left, v_x_bcs.right, v, Gv, 𝜈, u, v, w, T, S, t, iteration, grid)
    apply_bcs!(arch, Val(coord), Bx, By, Bz, w_x_bcs.left, w_x_bcs.right, w, Gw, 𝜈, u, v, w, T, S, t, iteration, grid)
    apply_bcs!(arch, Val(coord), Bx, By, Bz, T_x_bcs.left, T_x_bcs.right, T, GT, κ, u, v, w, T, S, t, iteration, grid)
    apply_bcs!(arch, Val(coord), Bx, By, Bz, S_x_bcs.left, S_x_bcs.right, S, GS, κ, u, v, w, T, S, t, iteration, grid)

    return nothing
end

# Do nothing if both boundary conditions are default.
apply_bcs!(::CPU, ::Val{:x}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing
apply_bcs!(::CPU, ::Val{:y}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing
apply_bcs!(::CPU, ::Val{:z}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing

apply_bcs!(::GPU, ::Val{:x}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing
apply_bcs!(::GPU, ::Val{:y}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing
apply_bcs!(::GPU, ::Val{:z}, Bx, By, Bz, left_bc::BC{<:Default, T}, right_bc::BC{<:Default, T}, args...) where {T} = nothing

# First, dispatch on coordinate.
apply_bcs!(arch, ::Val{:x}, Bx, By, Bz, args...) = @launch device(arch) apply_x_bcs!(args..., threads=(Tx, Ty), blocks=(Bx, By, Bz))
apply_bcs!(arch, ::Val{:y}, Bx, By, Bz, args...) = @launch device(arch) apply_y_bcs!(args..., threads=(Tx, Ty), blocks=(Bx, By, Bz))
apply_bcs!(arch, ::Val{:z}, Bx, By, Bz, args...) = @launch device(arch) apply_z_bcs!(args..., threads=(Tx, Ty), blocks=(Bx, By))

"Apply a top and/or bottom boundary condition to variable ϕ. Note that this kernel
MUST be launched with blocks=(Bx, By). If launched with blocks=(Bx, By, Bz), the
boundary condition will be applied Bz times!"
function apply_z_bcs!(top_bc, bottom_bc, ϕ, Gϕ, κ, u, v, w, T, S, t, iteration, grid::RegularCartesianGrid)
    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            apply_z_top_bc!(top_bc, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j)
            apply_z_bottom_bc!(bottom_bc, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j)
        end
    end
    @synchronize
end


text:
using TestFunctionRunner
TestFunctionRunner.@run(paths = ["../simplefail/SampleSimpleFail"])


text:
function bo_status_window!()
    win    = Gtk.Window("Gallery Status")
    vbox   = Gtk.Box(:v)
    txtbox = GtkReactive.textarea("Starting Bayesian Optimization")
    push!(win, vbox)
    push!(vbox, txtbox)
    Gtk.showall(win)
    dict   = Dict("window"=>win, "vbox"=>vbox)
    GtkReactive.gc_preserve(win, dict)
    dict
end

function bo_status_update!(dict::Dict, message::String)
    txtbox = GtkReactive.textarea(message)
    push!(dict["vbox"], txtbox)
end

function bo_status_destroy!(dict::Dict)
    GtkReactive.destroy(dict["window"])
end