text:
### A Pluto.jl notebook ###
# v0.11.10

using Markdown
using InteractiveUtils

# ╔═╡ 758de036-ecdd-11ea-3bb5-010d259ca80d
begin
	using HTTP, Pipe
end

# ╔═╡ 95053d12-ecdd-11ea-3cb2-27f60fc2eb5a
md"""
    downloader(url::String, dest::String, fname::String)

This function takes a url and a destination and downloads the file there with a name specified.\n
    url: The url to be downloaded
    dest: Path
    fname: File name
    ```julia
    downloader("http://www.julialang.org/", "/tmp/","index.html")
    ```
"""

# ╔═╡ 94eb9e8e-ecdd-11ea-345c-a121388fc870
function downloader(url::String, dest::String, fname::String)
    HTTP.open(:GET, url) do http
        open(dest*fname, "w") do file
            write(file, http)
        end
    end
    return dest*fname
end

# ╔═╡ 94d3f504-ecdd-11ea-33a6-0787b8e47ccd
md"""
An attempt to have an easy to download list of useful datasets in Machine Learning.
This list was mostly taken from the python fastai documentation and modified for Julia. [Source](https://course.fast.ai/datasets).
Some others were added from the Pytorch documentation. [Source](https://pytorch.org/docs/stable/torchvision/datasets.html)
Please refer to the sources for citations if I missed anything.
"""

# ╔═╡ 94bcd022-ecdd-11ea-35e2-9f09d7ea0cd6
dataset_list = Dict(

    # Computer Vision
        "mnist" => "https://s3.amazonaws.com/fast-ai-imageclas/mnist_png.tgz",
        "cifar10" => "https://s3.amazonaws.com/fast-ai-imageclas/cifar10.tgz",
        "cifar100" => "https://s3.amazonaws.com/fast-ai-imageclas/cifar100.tgz",
        "birds" => "https://s3.amazonaws.com/fast-ai-imageclas/CUB_200_2011.tgz",
        "caltech101" => "https://s3.amazonaws.com/fast-ai-imageclas/caltech_101.tar.gz",
        "pets" => "https://s3.amazonaws.com/fast-ai-imageclas/oxford-iiit-pet.tgz",
        "flowers" => "https://s3.amazonaws.com/fast-ai-imageclas/oxford-102-flowers.tgz",
        "food" => "https://s3.amazonaws.com/fast-ai-imageclas/food-101.tgz",
        "cars" => "https://s3.amazonaws.com/fast-ai-imageclas/stanford-cars.tgz",
        "imagenette" => "https://s3.amazonaws.com/fast-ai-imageclas/imagenette.tgz",
        "imagenette320" => "https://s3.amazonaws.com/fast-ai-imageclas/imagenette-320.tgz",
        "imagenette160" => "https://s3.amazonaws.com/fast-ai-imageclas/imagenette-160.tgz",
        "imagewoof" => "https://s3.amazonaws.com/fast-ai-imageclas/imagewoof.tgz",
        "imagewoof320" => "https://s3.amazonaws.com/fast-ai-imageclas/imagewoof-320.tgz",
        "imagewoof160" => "https://s3.amazonaws.com/fast-ai-imageclas/imagewoof-160.tgz",
    # NLP
        "imdb"=>"https://s3.amazonaws.com/fast-ai-nlp/imdb.tgz",
        "wikitext103"=>"https://s3.amazonaws.com/fast-ai-nlp/wikitext-103.tgz",
        "wikitext2"=>"https://s3.amazonaws.com/fast-ai-nlp/wikitext-2.tgz",
        "wmt"=>"https://s3.amazonaws.com/fast-ai-nlp/giga-fren.tgz",
        "ag"=>"https://s3.amazonaws.com/fast-ai-nlp/ag_news_csv.tgz",
        "amazon"=>"https://s3.amazonaws.com/fast-ai-nlp/amazon_review_full_csv.tgz",
        "p-amazon"=>"https://s3.amazonaws.com/fast-ai-nlp/amazon_review_polarity_csv.tgz",
        "dbpedia"=>"https://s3.amazonaws.com/fast-ai-nlp/dbpedia_csv.tgz",
        "sogou"=>"https://s3.amazonaws.com/fast-ai-nlp/sogou_news_csv.tgz",
        "yahoo"=>"https://s3.amazonaws.com/fast-ai-nlp/yahoo_answers_csv.tgz",
        "yelp"=>"https://s3.amazonaws.com/fast-ai-nlp/yelp_review_full_csv.tgz",
        "p-yelp"=>"https://s3.amazonaws.com/fast-ai-nlp/yelp_review_polarity_csv.tgz",
    # Image Localization
        "camvid"=>"https://s3.amazonaws.com/fast-ai-imagelocal/camvid.tgz",
        "pascal"=>"https://s3.amazonaws.com/fast-ai-imagelocal/pascal-voc.tgz",
    # Video
        "hmbd"=>"http://serre-lab.clps.brown.edu/wp-content/uploads/2013/10/hmdb51_org.rar", #H. Kuehne, H. Jhuang, E. Garrote, T. Poggio, and T. Serre. HMDB: A Large Video Database for Human Motion Recognition. ICCV, 2011
        "ucf"=>"https://www.crcv.ucf.edu/data/UCF101/UCF101.rar", #Khurram Soomro, Amir Roshan Zamir and Mubarak Shah, UCF101: A Dataset of 101 Human Action Classes From Videos in The Wild., CRCV-TR-12-01, November, 2012.
        "kinetics700" => "https://storage.googleapis.com/deepmind-media/Datasets/kinetics700.tar.gz", #https://arxiv.org/abs/1907.06987
        "kinetics600" => "https://storage.googleapis.com/deepmind-media/Datasets/kinetics600.tar.gz", #https://arxiv.org/abs/1808.01340
        "kinetics400" => "https://storage.googleapis.com/deepmind-media/Datasets/kinetics400.tar.gz" #https://arxiv.org/abs/1705.06950

)

# ╔═╡ 94a5d226-ecdd-11ea-3824-570792db6b99
md"""
Helper function to return a list of existing datasets
"""

# ╔═╡ 948a3c48-ecdd-11ea-335b-55571ecea719
alldatasets()= keys(dataset_list)

# ╔═╡ 1a473232-ecde-11ea-1b7e-630ba6c299e9
md"""
getext(st::String)

Helper function to get the extension of the datasets.
Takes into account patterns like .tar, .tar.gz , .zip etc..
"""

# ╔═╡ 1a3164d4-ecde-11ea-04e0-fd1c62efa36b
function getext(st::String)
	@pipe split(st, "/")[end] |> split(_, ".")[2:end] |> "."*join(_, ".")
end

# ╔═╡ 4ec7faae-ecee-11ea-3f38-03b3208a6a87
st = "https:Datasets/kinetics400.zip"

# ╔═╡ 5d0bbd76-ecee-11ea-3c9b-81ac883c4874
getext(st)

# ╔═╡ 2d9a6638-ecde-11ea-3c68-a187bba663c0
md"""
get_data(name::String, path::String)

Main function to download a standard dataset.
Takes input as name from the above list, the destination and a boolean to check if the user wants to extract the data.
The filename is inferred from the destination and the url.

```julia
get_data("mnist", "/tmp/", 1)
```

If not found, the list of available datasets is returned.
"""

# ╔═╡ 2d81ca74-ecde-11ea-3f54-db31ba1d8bf6
function get_data(name::String, path::String)

    if name in alldatasets()
        @info "Downloading $name"
        url = dataset_list[name]
        finfile = downloader(url, path, name*getext(url))
        # finfile = "/tmp/mnist.tgz"
        @info "Done downloading"
    else
        @info "Please choose something from here"; @info alldatasets()
    end
end

# ╔═╡ 2d6a7e82-ecde-11ea-0b1d-a3be01d621dc

# ╔═╡ 1a199d9a-ecde-11ea-0dc3-8510834a38b9

# ╔═╡ 948a184e-ecdd-11ea-1ba2-f9e84085759a

# ╔═╡ 94606696-ecdd-11ea-23ec-cd16b1e4eab3

# ╔═╡ Cell order:
# ╠═758de036-ecdd-11ea-3bb5-010d259ca80d
# ╟─95053d12-ecdd-11ea-3cb2-27f60fc2eb5a
# ╠═94eb9e8e-ecdd-11ea-345c-a121388fc870
# ╟─94d3f504-ecdd-11ea-33a6-0787b8e47ccd
# ╠═94bcd022-ecdd-11ea-35e2-9f09d7ea0cd6
# ╠═94a5d226-ecdd-11ea-3824-570792db6b99
# ╠═948a3c48-ecdd-11ea-335b-55571ecea719
# ╟─1a473232-ecde-11ea-1b7e-630ba6c299e9
# ╠═1a3164d4-ecde-11ea-04e0-fd1c62efa36b
# ╠═4ec7faae-ecee-11ea-3f38-03b3208a6a87
# ╠═5d0bbd76-ecee-11ea-3c9b-81ac883c4874
# ╟─2d9a6638-ecde-11ea-3c68-a187bba663c0
# ╠═2d81ca74-ecde-11ea-3f54-db31ba1d8bf6
# ╠═2d6a7e82-ecde-11ea-0b1d-a3be01d621dc
# ╠═1a199d9a-ecde-11ea-0dc3-8510834a38b9
# ╠═948a184e-ecdd-11ea-1ba2-f9e84085759a
# ╠═94606696-ecdd-11ea-23ec-cd16b1e4eab3


text:
"""
Capital Account

capital can be

"""

struct Capital :> crAccount

end

ext: jl

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using RobotDynamics
using StaticArrays
using Test
using BenchmarkTools
using ForwardDiff, FiniteDiff

using RobotDynamics: @autodiff, KnotPoint, getdata,
    StaticReturn, InPlace, ForwardAD, FiniteDifference
const RD = RobotDynamics

function test_model(model)
    n,m = size(model)
    xdot = zeros(n)
    x = @SVector randn(n)
    u = @SVector randn(m)
    t,dt = 0.0, 0.1
    zs = KnotPoint(x,u,t,dt)
    z = KnotPoint(Vector(x),Vector(u),t,dt)
    J = zeros(n, n + m)
    J0 = zeros(n, n + m)

    RD.evaluate!(model, xdot, z)
    @test xdot ≈ RD.evaluate(model, z)

    RD.evaluate!(model, xdot, x, u, RD.getparams(z))
    @test xdot ≈ RD.evaluate(model, x, u, RD.getparams(z))

    if model isa RD.ContinuousDynamics
        RobotDynamics.dynamics!(model, xdot, x, u, t)
        @test xdot ≈ RobotDynamics.dynamics(model, x, u, t)
    elseif model isa RD.DiscreteDynamics
        xn = xdot
        RobotDynamics.discrete_dynamics!(model, xn, x, u, t, dt)
        @test xn ≈ RobotDynamics.discrete_dynamics(model, x, u, t, dt)
    end

    RobotDynamics.jacobian!(StaticReturn(), RD.UserDefined(), model, J0, xdot, z)
    RobotDynamics.jacobian!(StaticReturn(), ForwardAD(), model, J, xdot, z)
    @test J ≈ J0
    RobotDynamics.jacobian!(InPlace(), RD.UserDefined(), model, J0, xdot, z)
    RobotDynamics.jacobian!(InPlace(), ForwardAD(), model, J, xdot, z)
    @test J ≈ J0
    RobotDynamics.jacobian!(StaticReturn(), FiniteDifference(), model, J, xdot, z)
    @test J ≈ J0 atol = 1e-6
    RobotDynamics.jacobian!(InPlace(), FiniteDifference(), model, J, xdot, z)
    @test J ≈ J0 atol = 1e-6

    allocs = 0
    if model isa RD.ContinuousDynamics
        allocs += @allocated RobotDynamics.dynamics!(model, xdot, x, u, t)
        allocs += @allocated RobotDynamics.dynamics(model, x, u, t)
    else
        allocs += @allocated RobotDynamics.discrete_dynamics!(model, xdot, x, u, t, dt)
        allocs += @allocated RobotDynamics.discrete_dynamics(model, x, u, t, dt)
    end
    allocs += @allocated RobotDynamics.evaluate!(model, xdot, z)
    allocs += @allocated RobotDynamics.evaluate(model, zs)
    allocs += @allocated RobotDynamics.evaluate!(model, xdot, x, u, RD.getparams(z))
    allocs += @allocated RobotDynamics.evaluate(model, x, u, RD.getparams(z))
    allocs += @allocated RobotDynamics.jacobian!(StaticReturn(), RD.UserDefined(), model, J, xdot, zs)
    allocs += @allocated RobotDynamics.jacobian!(InPlace(), RD.UserDefined(), model, J, xdot, zs)
    allocs += @allocated RobotDynamics.jacobian!(StaticReturn(), ForwardAD(), model, J, xdot, zs)
    allocs += @allocated RobotDynamics.jacobian!(InPlace(), ForwardAD(), model, J, xdot, z)
    allocs += @allocated RobotDynamics.jacobian!(StaticReturn(), FiniteDifference(), model, J, xdot, zs)
    allocs += @allocated RobotDynamics.jacobian!(InPlace(), FiniteDifference(), model, J, xdot, z)
    @test allocs == 0
end

function test_error_allocs(model)
    n,m = size(model)
    z1 = KnotPoint(randn(model)..., 0.0, 0.1)
    z2 = KnotPoint(randn(model)..., 0.1, 0.1)
    J1 = zeros(n, n+m)
    J2 = copy(J1)
    y1 = zeros(n)
    y2 = copy(y1)
    RD.dynamics_error!(model, y2, y1, z2, z1)
    @test y2 ≈ RD.dynamics_error(model, z2, z1)

    allocs = 0
    allocs += @allocated RD.dynamics_error_jacobian!(StaticReturn(), RD.UserDefined(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), RD.UserDefined(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(StaticReturn(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(StaticReturn(), FiniteDifference(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), FiniteDifference(), model, J2, J1, y2, y1, z2, z1)

    z1 = KnotPoint{n,m}(Vector(z1.z), z1.t, z1.dt)
    z2 = KnotPoint{n,m}(Vector(z2.z), z2.t, z2.dt)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), RD.UserDefined(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
    allocs += @allocated RD.dynamics_error_jacobian!(InPlace(), FiniteDifference(), model, J2, J1, y2, y1, z2, z1)
    @test allocs == 0
end

##
abstract type DI{D} <: RobotDynamics.ContinuousDynamics where {D} end
RobotDynamics.state_dim(::DI{D}) where {D} = 2D
RobotDynamics.control_dim(::DI{D}) where {D} = D

@generated function RobotDynamics.dynamics(fun::DI{D}, x, u, t) where {D}
    N, M = 2D, D
    vel = [:(x[$i]) for i = M+1:N]
    us = [:(u[$i]) for i = 1:M]
    :(SVector{$N}($(vel...), $(us...)))
end

@generated function RobotDynamics.dynamics!(fun::DI{D}, y, x, u, t) where {D}
    N, M = 2D, D
    vel = [:(y[$(i - M)] = x[$i]) for i = M+1:N]
    us = [:(y[$(i + M)] = u[$i]) for i = 1:M]
    quote
        $(Expr(:block, vel...))
        $(Expr(:block, us...))
        return nothing
    end
end

@generated function RobotDynamics.jacobian!(fun::DI{D}, J, y, x, u, t) where {D}
    N, M = 2D, D
    jac = [:(J[$i, $(i + M)] = 1.0) for i = 1:N]
    quote
        J .= 0
        $(Expr(:block, jac...))
        return nothing
    end
end


text:
abstract type KayRuleModel <: MixingRule end

struct KayRuleParam <: EoSParam
end

@newmodelsimple KayRule KayRuleModel KayRuleParam

export KayRule

function KayRule(components::Vector{String}; activity=nothing, userlocations::Vector{String}=String[], activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    #params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    #acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = KayRuleParam()
    model = KayRule(packagedparams, verbose=verbose)
    return model
end

function mixing_rule(model::ABCubicModel,V,T,z,mixing_model::KayRuleModel,α,a,b,c)
    n = sum(z)
    invn = (one(n)/n)
    invn2 = invn^2
    #b̄ = (dot(z,Symmetric(b.^(1/3)),z) * invn2).^3
    #ā = √(dot(z,Symmetric(a .* sqrt.(α*α') ./ b).^2,z)* invn2)  * b̄
    b̄ = zero(eltype(z))
    ab2 = zero(T+first(z))
    for i in 1:length(z)
        zi = z[i]
        zi2 = zi*zi
        αi = α[i]
        b̄ += cbrt(b[i,i])*zi2
        abij2 = (a[i,i]*αi/b[i,i])^2
        ab2 += abij2*zi2
        for j in 1:(i-1)
            zij =zi*z[j]
            abij2 = (a[i,j]*sqrt(αi*α[j])/b[i,j])^2
            ab2 += 2*abij2*zij
            b̄ += 2*cbrt(b[i,j])*zij
        end
    end
    b̄ = (b̄*invn2)^3
    ā = sqrt(ab2*invn2)*b̄
    c̄ = dot(z,c)*invn
    return ā,b̄,c̄
end

is_splittable(::KayRule) = false


text:
struct Assets

    GRID :: Grid
    Dx :: SparseMatrixCSC
    Dy :: SparseMatrixCSC
    Dxy :: SparseMatrixCSC
    Dxx :: SparseMatrixCSC
    Dyy :: SparseMatrixCSC
    Dxxx :: SparseMatrixCSC
    Dyyy :: SparseMatrixCSC
    Dxxy :: SparseMatrixCSC
    Dxyy :: SparseMatrixCSC
    Ds :: SparseMatrixCSC
    Dss :: SparseMatrixCSC
    DIFF_lnn :: LinearLeftHandSide
    DIFF_lnTe :: LinearLeftHandSide
    DIFF_lnTi :: LinearLeftHandSide
    DIFF_u :: LinearLeftHandSide
    DIFF_w :: LinearLeftHandSide
    DIFF_A :: LinearLeftHandSide
    HHOLTZ :: LinearLeftHandSide
    P0 :: SparseMatrixCSC
    P1 :: SparseMatrixCSC
    P2 :: SparseMatrixCSC
    P3 :: SparseMatrixCSC
    R1 :: SparseMatrixCSC
    R2 :: SparseMatrixCSC
    R3 :: SparseMatrixCSC
    LAM :: SparseMatrixCSC
    DCHLT1 :: SparseMatrixCSC
    DCHLT2 :: SparseMatrixCSC
    DCHLT3 :: SparseMatrixCSC
    NMANN1 :: SparseMatrixCSC
    NMANN2 :: SparseMatrixCSC
    NMANN3 :: SparseMatrixCSC
    FLXAVG :: SparseMatrixCSC
    TRGT :: Vector{Float64}
    Sn :: Vector{Float64}
    STe :: Vector{Float64}
    STi :: Vector{Float64}
    params :: Vector{Float64}
    dt :: Float64
    N_subcycle :: Int64

end

function Assets(geoinit::Function, phys_ref::Vector{Float64}, Nz::Int64, ds::Float64, dt::Float64; k=0.2, w=0.66, mx=5, my=5, N_subcycle=100)

    # return geoinit()
    grd, bdry, lcfs_pts, trgt_splitter, srcs, bx, by, bz = geoinit()
    params = dimensionless_parameters(phys_ref...)

    P0, P1, P2, P3 = bdry[1]
    R1, R2, R3 = bdry[2]
    DCHLT1, DCHLT2, DCHLT3 = bdry[3]
    NMANN1, NMANN2, NMANN3 = bdry[4]
    Sn, STe, STi = srcs

    # generate stencils
    Mxinv = stencil1d(mx)
    Myinv = stencil1d(my)
    MxyinvT = stencil2d(mx, my)

    # generate derivatives
    Dx = derivative_matrix(1, 0, Mxinv, Myinv, grd)
    Dy = derivative_matrix(0, 1, Mxinv, Myinv, grd)
    Dxy = derivative_matrix(1, 1, Mxinv, Myinv, grd)
    Dxx = derivative_matrix(2, 0, Mxinv, Myinv, grd)
    Dyy = derivative_matrix(0, 2, Mxinv, Myinv, grd)
    Dxxx = derivative_matrix(3, 0, Mxinv, Myinv, grd)
    Dyyy = derivative_matrix(0, 3, Mxinv, Myinv, grd)
    Dxxy = derivative_matrix(2, 1, Mxinv, Myinv, grd)
    Dxyy = derivative_matrix(1, 2, Mxinv, Myinv, grd)
    Ds, Dss = fieldline_derivatives(bx, by, bz, ds, mx, my, MxyinvT, Nz, grd)

    @info "Generated derivative matrices."

    # compute LAM
    LAM = Diagonal(1 ./ (1 .+ dt*diag(P1 + P2 + P3)))

    # setup linear problems
    de2 = params[8]
    Ldiff = I - k*(Dxx + Dyy)
    Lhelm = I - de2*(Dxx + Dyy)
    DIFF_lnn = LinearLeftHandSide(P0*Ldiff + NMANN1 + NMANN2 + NMANN3, w)
    DIFF_lnTe = LinearLeftHandSide(P0*Ldiff + NMANN1 + NMANN2 + NMANN3, w)
    DIFF_lnTi = LinearLeftHandSide(P0*Ldiff + NMANN1 + NMANN2 + NMANN3, w)
    DIFF_u = LinearLeftHandSide(P0*Ldiff + NMANN1 + NMANN2 + DCHLT3, w)
    DIFF_w = LinearLeftHandSide(P0*Ldiff + DCHLT1 + DCHLT2 + DCHLT3, w)
    DIFF_A = LinearLeftHandSide(P0*Ldiff + DCHLT1 + DCHLT2 + DCHLT3, w)
    HHOLTZ = LinearLeftHandSide(P0*Lhelm + DCHLT1 + DCHLT2 + DCHLT3, w)

    @info "Setup linear problems."

    # compute vector for averaging over the LCFS
    rows = SparseMatrixCSC[]
    for i=1:size(lcfs_pts)[2]
        dat, js = interpolation_row(lcfs_pts[:,i]..., mx, my, MxyinvT, grd)
        is = ones(Int64, length(js))
        row = sparse(is, js, dat, 1, grd._Nx*grd._Ny) * transpose(grd.Proj)
        append!(rows, [row])
    end
    FLXAVG = +(rows...) / size(lcfs_pts)[2]

    # define TRGT
    TRGT = f_to_grid((x, y) -> sign(poly_dist(x, y, trgt_splitter)[1]), grd)

    return Assets(grd, Dx, Dy, Dxy, Dxx, Dyy, Dxxx, Dyyy, Dxxy, Dxyy, Ds, Dss,
                     DIFF_lnn, DIFF_lnTe, DIFF_lnTi, DIFF_u, DIFF_w, DIFF_A,
                     HHOLTZ, P0, P1, P2, P3, R1, R2, R3, LAM, DCHLT1, DCHLT2,
                     DCHLT3, NMANN1, NMANN2, NMANN3, FLXAVG, TRGT, Sn, STe,
                     STi, params, dt, N_subcycle)
end

function simple_solovev_assets(Nx, Ny, A, ep, ka, del, L, R0, n0, T0, B0)

    phys_ref = [ep*R0, R0, n0, T0, B0]

    Assets(phys_ref, 1, 0.05, 0.001) do

        psi, bx, by, bz = solovev_flux_function(A, del, ep, ka, 100)
        bxp(x,y) = bx(x,y) / norm([bx(x,y), by(x,y)])
        byp(x,y) = by(x,y) / norm([bx(x,y), by(x,y)])
        psi_in = psi(1 - ep + L, 0)
        psi_out = 0

        inner_flux_surface = Barrier() do
            rmap(x, y) = trace_reflection(x, y, psi, bxp, byp, psi_in, 0.005)
            psi, psi_in, rmap
        end

        outer_flux_surface = Barrier() do
            rmap(x, y) = trace_reflection(x, y, psi, bxp, byp, psi_out, 0.001)
            psi, psi_out, rmap
        end

        bars = [inner_flux_surface, outer_flux_surface, Everywhere()]
        bigdeltas = [1, -1, 1] * abs(psi_in) * 1.5
        grd = Grid(Nx, Ny, 1) do
            r0 = [0.01, -(ka*ep + 5*L)]
            r1 = [1 + ep+L, ka*ep + 5*L]
            (x,y) -> check_if_inside(x, y, bigdeltas, bars), r0, r1
        end

        # return grd, bars, psi

        @info "Generated Grid"

        smalldeltas = bigdeltas / 2.0
        qs = ones(3) * 0.01
        bdry = boundary_operators(3, 3, stencil2d(3, 3), smalldeltas, bars, qs, grd)

        @info "Generated boundary operators."

        lcfs_pts = zeros(2, 30)
        ts = LinRange(0, 2*pi, 31)
        ts = ts[1:end-1]
        alpha = asin(del)
        for i=1:30
            lcfs_pts[1,i] = 1 + (ep-L/2)*cos(ts[i] + alpha*sin(ts[i]))
            lcfs_pts[2,i] = (ep-L/2)*ka*sin(ts[i])
        end

        # trgt_splitter doesn't matter because there is no limiter or diverter here
        trgt_splitter = Float64[0 1; 0 1]

        # define constant source for now
        Sn = f_to_grid((x,y) -> 1/(100*0.001), grd)
        STe = f_to_grid((x,y) -> 1/(100*0.001), grd)
        STi = f_to_grid((x,y) -> 1/(100*0.001), grd)

        grd, bdry, lcfs_pts, trgt_splitter, [Sn, STe, STi], bx, by, bz
    end

end


text:
@test collect(partitions(4)) ==  Any[[4], [3,1], [2,2], [2,1,1], [1,1,1,1]]
@test collect(partitions(8,3)) == Any[[6,1,1], [5,2,1], [4,3,1], [4,2,2], [3,3,2]]
@test collect(partitions(8, 1)) == Any[[8]]
@test collect(partitions(8, 9)) == []
@test collect(partitions([1,2,3])) == Any[Any[[1,2,3]], Any[[1,2],[3]], Any[[1,3],[2]], Any[[1],[2,3]], Any[[1],[2],[3]]]
@test collect(partitions([1,2,3,4],3)) == Any[Any[[1,2],[3],[4]], Any[[1,3],[2],[4]], Any[[1],[2,3],[4]],
                                              Any[[1,4],[2],[3]], Any[[1],[2,4],[3]], Any[[1],[2],[3,4]]]
@test collect(partitions([1,2,3,4],1)) == Any[Any[[1, 2, 3, 4]]]
@test collect(partitions([1,2,3,4],5)) == []

@inferred first(partitions(4))
@inferred first(partitions(8,3))
@inferred first(partitions([1,2,3]))
@inferred first(partitions([1,2,3,4],3))

@test isa(collect(partitions(4)), Vector{Vector{Int}})
@test isa(collect(partitions(8,3)), Vector{Vector{Int}})
@test isa(collect(partitions([1,2,3])), Vector{Vector{Vector{Int}}})
@test isa(collect(partitions([1,2,3,4], 3)), Vector{Vector{Vector{Int}}})

@test length(partitions(0)) == 1
@test length(partitions(-1)) == 0
@test length(collect(partitions(30))) == length(partitions(30))
@test length(collect(partitions(90,4))) == length(partitions(90,4))
@test length(collect(partitions('a':'h'))) == length(partitions('a':'h'))
@test length(collect(partitions('a':'h',5))) == length(partitions('a':'h',5))

# integer_partitions
@test integer_partitions(0) == []
@test integer_partitions(5) == Any[[1, 1, 1, 1, 1], [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]]
@test_throws DomainError integer_partitions(-1)

@test_throws ArgumentError prevprod([2,3,5],Int128(typemax(Int))+1)
@test prevprod([2,3,5],30) == 30
@test prevprod([2,3,5],33) == 32

# noncrossing partitions
let nc4 = ncpartitions(4)
    @test nc4 == Any[Any[[1],[2],[3],[4]], Any[[1],[2],[3,4]], Any[[1],[2,3],[4]], Any[[1],[2,4],[3]], Any[[1],[2,3,4]],
                     Any[[1,2],[3],[4]], Any[[1,2],[3,4]],
                     Any[[1,3],[2],[4]], Any[[1,4],[2],[3]], Any[[1,4],[2,3]],
                     Any[[1,2,3],[4]], Any[[1,3,4],[2]], Any[[1,2,4],[3]],
                     Any[[1,2,3,4]]]
    @test length(nc4) == catalannum(4)
end


text:
using BenchmarkTools: leaves
using PkgBenchmark: benchmarkgroup
using TransducersBenchmarksReports: loadresult
using DataFrames
using Statistics
using VegaLite

fanout(; kwargs...) =
    x -> NamedTuple{keys(kwargs)}(Ref(x) .|> Tuple(values(kwargs)))

macro fanout(names...)
    kwargs = [:($n=$n) for n in names]
    esc(:($fanout(; $(kwargs...))))
end

group = benchmarkgroup(loadresult())["gemm"]
data = map(leaves(group["fusedmul"])) do (key, trial)
    stats = @fanout(mean, median, minimum)(trial.times)
    impl, m = key
    return (impl=impl, m=m, stats...)
end

df = melt(DataFrame(data), [:impl, :m])
df[:m] = parse.(Int, df[:m])
sort!(df, :m)

dfnorm = join(
    df[df[:impl] .!= "blas", :],
    deletecols!(df[df[:impl] .== "blas", :], :impl),
    on = [:m, :variable],
    kind = :left,
    makeunique = true,
)
dfnorm[:normalized] = dfnorm[:value] ./ dfnorm[:value_1]

plt = dfnorm[dfnorm[:variable] .== :median, :] |> @vlplot(
    mark = {
        typ = :line,
        point = true,
    },
    x = {
        field = :m,
        typ = :quantitative,
        scale = {typ = :log},
    },
    y = {
        field = :normalized,
        typ = :quantitative,
        # scale = {typ = :log},
        axis = {title = "transducers / blas [median time]"},
    },
)


text:
using DFControl, Test
configure_pseudoset("test", joinpath(testdir, "testassets", "pseudos"), "localhost_test")
add_environment(DFControl.Client.environment_from_jobscript(joinpath(testdir, "testassets", "reference_job", "job.tt"), "localhost_test"), "test_default", "localhost_test")
# @test DFControl.Service.pseudos("test", "UPF")[:Ni] ==
#       Pseudo("Ni.pbesol-n-kjpaw_psl.0.1.UPF", joinpath(testdir, "testassets", "pseudos"), 41., 236.)