text:
"""
    mask(ctx, A)

"""
function mask(ctx, V)
    (μ⁰, _) = getproperty(ctx, :μ)
    μ⁰ * V
end

function mask(ctx, A::Union{TupleN{T},AbstractVector{T}}) where {T<:AbstractVector}
    (_, μ¹) = getproperty(ctx, :μ)
    μ¹ .* A
end

"""
    permanent(ctx, Φ, Ψ)

!!! note "Factor 2!"

"""
function permanent(ctx, Φ, Ψ)
    (τ⁻, τ⁺) = getproperty(ctx, :τ)
    map(getproperty(ctx, :τ)) do iter
        map(iter) do el
            ((el * Φ) .* Ψ) .+ (Φ .* (el * Ψ))
        end
    end
end

"""
    gradient(ctx, A, V, T, D)

Compute the gradient of T with Dirichlet boundary condition D.

"""
function gradient(::Dirichlet, ctx, A, V, T, D)
    ((δ⁻, δ⁺), (σ⁻, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(eachindex(A)) do i
        ((δ⁻[i] * ((σ⁺[i] * A[i]) .* T)) .- (σ⁻[i] * ((δ⁺[i] * A[i]) .* D))) ./ (σ⁻[i] * V)
    end
end

"""
    gradient(ctx, A, V, H, T, N)

Compute the gradient of T with Neumann boundary condition N.

!!! note

    Would it make sense to let N be a cell-centered vector-valued field?

"""
function gradient(::Neumann, ctx, A, V, H, T, N)
    ((δ⁻, δ⁺), (σ⁻, _)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(eachindex(A)) do i
        (2A[i] .* (δ⁻[i] * T) - (σ⁻[i] * ((δ⁺[i] * A[i]) .* H .* N))) ./ (σ⁻[i] * V)
    end
end

"""
    gradient(ctx, A, V, H, T, N)

Compute (twice) the volume-weighted gradient of T with Neumann boundary condition N.

"""
function gradient(::Neumann, ctx, A, H, T, N)
    ((δ⁻, δ⁺), (σ⁻, _)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(eachindex(A)) do i
        (2A[i] .* (δ⁻[i] * T) - (σ⁻[i] * ((δ⁺[i] * A[i]) .* H .* N)))
    end
end

"""
    gradient(ctx, bc, A, V, T, D)

Compute the gradient of T.

"""
#function gradient(ctx, bc, A, V, T, D)
function gradient((; info)::Mixed, ctx, A, V, H, T, D, N)
    ((δ⁻, δ⁺), (σ⁻, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(eachindex(A)) do i
        ((δ⁻[i] * ((σ⁺[i] * A[i]) .* T)) .-
         (σ⁻[i] * ((δ⁺[i] * A[i]) .* ((info[i] .* D) .+ ((!).(info[i]) .* (T .+ (H .* N))))))) ./
        (σ⁻[i] * V)
    end
end

"""
    divergence(ctx, A, U)

Compute the volume-weighted divergence of U with Dirichlet boundary conditions.

!!! note "Factor 2!"

"""
function divergence(::Dirichlet, ctx, A, U)
    ((_, δ⁺), (_, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    mapreduce(+, eachindex(U)) do i
        (σ⁺[i] * A[i]) .* (δ⁺[i] * U[i])
    end
end

"""
    divergence(bc, ctx, A, U, W)

Compute the volume-weighted divergence of U with Neumann boundary conditions.

!!! note "Factor 2!"

!!! warning

    It would make a lot of sense to let W be a cell-centered vector-valued field.

"""
function divergence(::Neumann, ctx, A, U, W)
    (_, δ⁺) = getproperty(ctx, :δ)
    mapreduce(+, eachindex(U)) do i
        (2δ⁺[i] * (A[i] .* U[i])) .- ((2δ⁺[i] * A[i]) .* W)
    end
end

"""
    divergence(ctx, bc, A, U, W)

!!! note "Factor 2!"

"""
function divergence((; info)::Mixed, ctx, A, U, W)
    ((_, δ⁺), (_, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    mapreduce(+, eachindex(A)) do i
        (2δ⁺[i] * (A[i] .* U[i])) .-
        ((δ⁺[i] * A[i]) .* ((info[i] .* (σ⁺[i] * U[i])) .+ ((!).(info[i]) .* 2W)))
    end
end

"""
    strainrate(ctx, A, V, U, W)

!!! note

    For now, only diagonal terms (and extra factor 2).
    Also, Dirichlet.

"""
function strainrate(::Dirichlet, ctx, A, V, U, W)
    ((δ⁻, δ⁺), (σ⁻, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(product(eachindex(A), eachindex(A))) do (i, j)
        if i == j
            ((2δ⁺[i] * (A[i] .* U[i])) .- ((δ⁺[i] * A[i]) .* (σ⁺[i] * W[i]))) ./ V
        else
            (2δ⁻[j] * ((σ⁻[i] * (σ⁺[j] * A[j])) .* U[i])
             .- (σ⁻[j] * ((σ⁻[i] * (δ⁺[j] * A[j])) .* W[i]))) ./ (σ⁻[j] * (σ⁻[i] * V)) .+
            (2δ⁻[i] * ((σ⁻[j] * (σ⁺[i] * A[i])) .* U[j])
             .- (σ⁻[i] * ((σ⁻[j] * (δ⁺[i] * A[i])) .* W[j]))) ./ (σ⁻[i] * (σ⁻[j] * V))
        end
    end
end

"""
    divergence(ctx, A, E)

!!! note

    Also, Dirichlet.

"""
function divergence(::Dirichlet, ctx, A, E::AbstractMatrix)
    ((δ⁻, δ⁺), (σ⁻, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    map(eachindex(A)) do i
        mapreduce(+, eachindex(A)) do j
            if i == j
                4A[i] .* (δ⁻[i] * E[i, i])
            else
                (σ⁻[i] * (σ⁺[j] * A[j])) .* (δ⁺[j] * E[i, j])
            end
        end
    end
end

function dissipation(ctx, A, V, T)
    ((δ⁻, δ⁺), (σ⁻, σ⁺)) = getproperty.(Ref(ctx), (:δ, :σ))
    mapreduce(+, eachindex(A)) do i
        -(σ⁺[i] * ((δ⁻[i] * ((σ⁺[i] * A[i]) .* T)) .^ 2 ./ (σ⁻[i] * V)))
    end
end


text:
using DataDeps

register(DataDep(
  "geographic-origin-music",
  "https://archive.ics.uci.edu/ml/datasets/Geographical+Original+of+Music",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00315/Geographical%20Original%20of%20Music.zip",
  "37bd14730ef7e4786e094421982ad32536298a9d84d875cd35d27423f14b409a",
  post_fetch_method=(path -> begin
    unpack(path)
    new_path = joinpath(dirname(path), "Geographical Original of Music/default_plus_chromatic_features_1059_tracks.txt")
    mv(new_path, path)

    UCIData.process_dataset(path,
        target_index=69,
        feature_indices=1:68,
    )
  end),
))

ext: jl

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
#=
# Plotting Graphs with `GraphMakie.jl`
## The `graphplot` Command
Plotting your first `AbstractGraph` from [`Graphs.jl`](https://github.com/JuliaGraphs/Graphs.jl)
is as simple as
=#
using CairoMakie
CairoMakie.activate!(type="png") # hide
set_theme!(resolution=(800, 400)) #hide
CairoMakie.inline!(true) # hide
using GraphMakie
using Graphs
import Random; Random.seed!(2) # hide

g = wheel_graph(10)
f, ax, p = graphplot(g)
hidedecorations!(ax); hidespines!(ax)
ax.aspect = DataAspect()
f # hide

#=
The `graphplot` command is a recipe which wraps several steps
- layout the graph in space using a layout function,
- create a `scatter` plot for the nodes and
- create a `linesegments` plot for the edges.

The default layout is `Spring()` from
[`NetworkLayout.jl`](https://github.com/JuliaGraphs/NetworkLayout.jl). The
layout attribute can be any function which takes an `AbstractGraph` and returns
a list of `Point{dim,Ptype}` (see [`GeometryBasics.jl`](https://github.com/JuliaGeometry/GeometryBasics.jl)
objects where `dim` determines the dimensionality of the plot.

Besides that there are some common attributes which are forwarded to the
underlying plot commands. See [`graphplot`](@ref).
=#
using GraphMakie.NetworkLayout

g = SimpleGraph(5)
add_edge!(g, 1, 2); add_edge!(g, 2, 4);
add_edge!(g, 4, 3); add_edge!(g, 3, 2);
add_edge!(g, 2, 5); add_edge!(g, 5, 4);
add_edge!(g, 4, 1); add_edge!(g, 1, 5);

## define some edge colors
edgecolors = [:black for i in 1:ne(g)]
edgecolors[4] = edgecolors[7] = :red

f, ax, p = graphplot(g, layout=Shell(),
                     node_color=[:black, :red, :red, :red, :black],
                     edge_color=edgecolors)

hidedecorations!(ax); hidespines!(ax)
ax.aspect = DataAspect()
f #hide

#=
We can interactively change the attributes as usual with Makie.
=#

fixed_layout(_) = [(0,0), (0,1), (0.5, 1.5), (1,1), (1,0)]
## set new layout
p.layout = fixed_layout; autolimits!(ax)
## change edge width & color
p.edge_width = 5.0
p.edge_color[][3] = :green;
p.edge_color = p.edge_color[] # trigger observable
f #hide

#=
## Adding Node Labels
=#
Random.seed!(2)
g = wheel_graph(10)

colors = [:black for i in 1:nv(g)]
colors[1] = :red

f, ax, p = graphplot(g,
                     nlabels=repr.(1:nv(g)),
                     nlabels_color=colors,
                     nlabels_align=(:center,:center))
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

# This is not very nice, lets change the offsets based on the `node_positions`

offsets = 0.15 * (p[:node_pos][] .- p[:node_pos][][1])
offsets[1] = Point2f(0, 0.3)
p.nlabels_offset[] = offsets
autolimits!(ax)
f # hide

# ## Adding Edge Labels
Random.seed!(42)
g = barabasi_albert(6, 2)

labels =  repr.(1:ne(g))

f, ax, p = graphplot(g, elabels=labels,
                     elabels_color=[:black for i in 1:ne(g)],
                     edge_color=[:black for i in 1:ne(g)])
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
The position of the edge labels is determined by several plot arguments.
All possible arguments are described in the docs of the [`graphplot`](@ref) function.

Basicially, each label is placed in the middle of the edge and rotated to match the edge rotation.
The rotaion for each label can be overwritten with the `elabels_rotation` argument.
=#
p.elabels_rotation[] = Vector{Union{Nothing, Float64}}(nothing, ne(g))
p.elabels_rotation[][5] = 0.0 # set absolute rotation angle for label 5
p.elabels_rotation[] = p.elabels_rotation[]
nothing #hide

#=
One can shift the label along the edge with the `elabels_shift` argument and determine the distance
in pixels using the `elabels_distance` argument.
=#
p.elabels_opposite[] = [1,2,8,6]

p.elabels_offset[] = [Point2f(0.0, 0.0) for i in 1:ne(g)]
p.elabels_offset[][5] = Point2f(-0.4,0)
p.elabels_offset[] = p.elabels_offset[]

p.elabels_shift[] = [0.5 for i in 1:ne(g)]
p.elabels_shift[][1] = 0.6
p.elabels_shift[][7] = 0.4
p.elabels_shift[] = p.elabels_shift[]

p.elabels_distance[] = zeros(ne(g))
p.elabels_distance[][8] = 15
p.elabels_distance[] = p.elabels_distance[]

f # hide

#=
## Indicate Edge Direction

It is possible to put arrows on the edges using the `arrow_show` parameter. This parameter
is `true` for `SimpleDiGraph` by default. The position and size of each arrowhead can be
change using the `arrow_shift` and `arrow_size` parameters.
=#
g = wheel_digraph(10)
arrow_size = [10+i for i in 1:ne(g)]
arrow_shift = range(0.1, 0.8, length=ne(g))
f, ax, p = graphplot(g; arrow_size, arrow_shift)
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
## Self edges

A self edge in a graph will be displayed as a loop.

!!! note
    Selfe edges are not possible in 3D plots yet.
=#
g = complete_graph(3)
add_edge!(g, 1, 1)
add_edge!(g, 2, 2)
add_edge!(g, 3, 3)
f, ax, p = graphplot(g)

hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

# It is possible to change the appearance using the `selfedge_` attributes:
p.selfedge_size = Dict(1=>Makie.automatic, 4=>3.6, 6=>0.5) #idx as in edges(g)
p.selfedge_direction = Point2f(0.3, 1)
p.selfedge_width = Any[Makie.automatic for i in 1:ne(g)]
p.selfedge_width[][4] = 0.6*π; notify(p.selfedge_width)
autolimits!(ax)
f # hide

#=
## Curvy edges

Curvy edges are possible using the low level interface of passing tangent
vectors and a `tfactor`. The tangent vectors can be `nothing` (straight line) or
two vectors per edge (one for src vertex, one for dst vertex). The `tfactor`
scales the distance of the bezier control point relative to the distance of src
and dst nodes. For real world usage see the [AST of a Julia function](@ref) example.
=#
using GraphMakie: plot_controlpoints!
g = complete_graph(3)
tangents = Dict(1 => ((1,1),(0,-1)),
                2 => ((0,1),(0,-1)),
                3 => ((0,-1),(1,0)))
tfactor = [0.5, 0.75, (0.5, 0.25)]
f, ax, p = graphplot(g; layout=SquareGrid(cols=3), tangents, tfactor,
                     arrow_size=20, arrow_show=true, edge_color=[:red, :green, :blue],
                     elabels="Edge ".*repr.(1:ne(g)), elabels_distance=10)
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
plot_controlpoints!(ax, p) # show control points for demonstration
f # hide

#=
## Edge waypoints
It is possible to specify waypoints per edge which needs to be crossed. See the
[Dependency Graph of a Package](@ref) example.

If the attribute `waypoint_radius` is `nothing` or `:spline` the waypoints will be crossed
using natural cubic spline interpolation. If the supply a radius the waypoints won't be reached,
instead they will be connected with straight lines which bend in the given radius around the
waypoints.
=#
set_theme!(resolution=(800, 800)) #hide
g = SimpleGraph(8); add_edge!(g, 1, 2); add_edge!(g, 3, 4); add_edge!(g, 5, 6); add_edge!(g, 7, 8)

waypoints = Dict(1 => [(.25,  0.25), (.75, -0.25)],
                 2 => [(.25, -0.25), (.75, -0.75)],
                 3 => [(.25, -0.75), (.75, -1.25)],
                 4 => [(.25, -1.25), (.75, -1.75)])
waypoint_radius = Dict(1 => nothing,
                       2 => 0,
                       3 => 0.05,
                       4 => 0.15)

f = Figure(); f[1,1] = ax = Axis(f)
using Makie.Colors # hide
for i in 3:4 #hide
    poly!(ax, Circle(Point2f(waypoints[i][1]), waypoint_radius[i]), color=RGBA(0.0,0.44705883,0.69803923,0.2)) #hide
    poly!(ax, Circle(Point2f(waypoints[i][2]), waypoint_radius[i]), color=RGBA(0.0,0.44705883,0.69803923,0.2)) #hide
end #hide

p = graphplot!(ax, g; layout=SquareGrid(cols=2, dy=-0.5),
               waypoints, waypoint_radius,
               nlabels=["","r = nothing (equals :spline)",
                        "","r = 0 (straight lines)",
                        "","r = 0.05 (in data space)",
                        "","r = 0.1"],
               nlabels_distance=30, nlabels_align=(:left,:center))

for i in 1:4 #hide
    scatter!(ax, waypoints[i], color=RGBA(0.0,0.44705883,0.69803923,1.0)) #hide
end #hide
xlims!(ax, (-0.1, 2.25)), hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
## Plot Graphs in 3D
If the layout returns points in 3 dimensions, the plot will be in 3D. However this is a bit
experimental. Feel free to file an issue if there are any problems.
=#
set_theme!(resolution=(800, 800)) #hide
g = smallgraph(:cubical)
elabels_shift = [0.5 for i in 1:ne(g)]
elabels_shift[[2,7,8,9]] .= 0.3
elabels_shift[10] = 0.25
graphplot(g; layout=Spring(dim=3, seed=5),
          elabels="Edge ".*repr.(1:ne(g)),
          elabels_textsize=12,
          elabels_opposite=[3,5,7,8,12],
          elabels_shift,
          elabels_distance=3,
          arrow_show=true,
          arrow_shift=0.9,
          arrow_size=15)

#=
Using [`JSServe.jl`](https://github.com/SimonDanisch/JSServe.jl) and [`WGLMakie.jl`](https://github.com/JuliaPlots/WGLMakie.jl)
we can also add some interactivity:
=#
using JSServe
Page(exportable=true, offline=true)
#
using WGLMakie
WGLMakie.activate!()
set_theme!(resolution=(800, 600))
g = smallgraph(:dodecahedral)
graphplot(g, layout=Spring(dim=3), node_size=100)


text:
function generate_batch_data(
    mountaincar::MountainCar,
    params::MountainCarParameters,
    num_episodes::Int64,
    horizon::Int64;
    policy = nothing,
    rng = nothing,
)::Array{MountainCarContTransition}
    data = []
    for episode = 1:num_episodes
        data = vcat(data, simulate_episode(mountaincar, params, horizon; policy=policy, rng = rng))
    end
    data = vcat(data, simulate_episode(mountaincar, params, horizon; policy = good_policy(mountaincar), rng = rng))
    data
end

function simulate_episode(
    mountaincar::MountainCar,
    params::MountainCarParameters,
    horizon::Int64;
    policy = nothing,
    rng = nothing,
)::Array{MountainCarContTransition}
    episode_data = []
    cont_state = init(mountaincar; cont = true)
    if isnothing(policy)
        policy = random_policy(mountaincar; rng = rng)
        cont_state = init(mountaincar; rng = rng, cont = true)
    end
    actions = getActions(mountaincar)
    for t = 1:horizon
        action = actions[policy[cont_state_to_idx(mountaincar, cont_state)]]
        cont_state_next, cost = step(mountaincar, cont_state, action, params; rng = rng)
        push!(
            episode_data,
            MountainCarContTransition(cont_state, action, cost, cont_state_next),
        )
        if checkGoal(mountaincar, cont_state_next)
            break
        end
        cont_state = cont_state_next
    end
    episode_data
end

function random_policy(mountaincar::MountainCar; rng = nothing)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization
    n_actions = 2
    if isnothing(rng)
        rand(1:n_actions, n_states)
    else
        rand(rng, 1:n_actions, n_states)
    end
end

function good_policy(mountaincar::MountainCar)
    value_iteration(mountaincar, vec(true_params))[1]
end

function preprocess_data(data::Array{MountainCarContTransition})
    n_actions = 2
    x_array::Array{Array{Array{Float64}}} = []
    x_next_array::Array{Array{Array{Float64}}} = []
    disp_array::Array{Array{Array{Float64}}} = []
    cost_array::Array{Array{Float64}} = []
    for a = 1:n_actions
        push!(x_array, [])
        push!(x_next_array, [])
        push!(disp_array, [])
        push!(cost_array, [])
    end

    for i = 1:length(data)
        transition = data[i]
        a = transition.action.id + 1
        push!(x_array[a], vec(transition.initial_state))
        push!(x_next_array[a], vec(transition.final_state))
        push!(disp_array[a], vec(transition.final_state) - vec(transition.initial_state))
        push!(cost_array[a], transition.cost)
    end
    x_array_matrices::Array{Matrix{Float64}} = [permutedims(hcat(x_subarray...)) for x_subarray in x_array]
    x_array_matrices, x_array, x_next_array, disp_array, cost_array
end

function eval_params(
    mountaincar::MountainCar,
    model::MountainCar,
    params::Array{Float64},
    horizon::Int64;
    num_eval = 10,
    rng = nothing,
)
    policy, _, _ = value_iteration(model, params)
    num_steps = 0.0
    for _ in 1:num_eval
        num_steps += length(
            simulate_episode(
                mountaincar,
                true_params,
                horizon;
                policy = policy,
                rng = rng,
            )
        )
    end
    num_steps/num_eval
end

function eval_params_bellman(
    mountaincar::MountainCar,
    model::MountainCar,
    params::Array{Float64},
    horizon::Int64;
    num_eval = 10,
    rng = nothing,
)
    policy, values, _ = value_iteration(model, params)
    num_steps = 0.0
    actions = getActions(mountaincar)
    for _ in 1:num_eval
        episode = simulate_episode(
            mountaincar,
            true_params,
            horizon;
            policy = policy,
            rng = rng,
        )
        model_steps = length(simulate_episode(
            model,
            unvec_params(params),
            horizon;
            policy = policy,
            rng = rng,
        ))
        println(model_steps, " ", values[cont_state_to_idx(model, init(mountaincar; cont = true))])
        for transition in episode
            successor = cont_state_to_idx(model, transition.final_state)
            action = actions[policy[cont_state_to_idx(model, transition.initial_state)]]
            predicted_successor = cont_state_to_idx(
                model, step(model, transition.initial_state, action, params)[1]
            )
            model_steps += values[successor] - values[predicted_successor]
        end
        model_steps += -values[cont_state_to_idx(model, episode[end].final_state)]
        num_steps += model_steps
    end
    num_steps/num_eval
end


text:
module TestTunedModels

using Distributed

using Test
using MLJTuning
using MLJBase
import ComputationalResources: CPU1, CPUProcesses, CPUThreads
using Random
Random.seed!(1234)
@everywhere using ..Models
using ..TestUtilities

N = 30
x1 = rand(N);
x2 = rand(N);
x3 = rand(N);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.4*rand(N);

m(K) = KNNRegressor(K=K)
r = [m(K) for K in 2:13]

@testset "constructor" begin
    @test_throws ErrorException TunedModel(model=first(r), tuning=Explicit(),
                                           measure=rms)
    @test_throws ErrorException TunedModel(tuning=Explicit(),
                                           range=r, measure=rms)
    @test_throws ErrorException TunedModel(model=42, tuning=Explicit(),
                                           range=r, measure=rms)
    @test_logs((:info, r"No measure"),
               TunedModel(model=first(r), tuning=Explicit(), range=r))
end

results = [(evaluate(model, X, y,
                     resampling=CV(nfolds=2),
                     measure=rms,
                     verbosity=0,)).measurement[1] for model in r]

@testset "measure compatibility check" begin
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=cross_entropy)
    @test_throws ArgumentError fit(tm, 0, X, y)
end

@testset_accelerated "basic fit" accel (exclude=[CPUThreads],) begin
    best_index = argmin(results)
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=[rms, l1], acceleration=accel)
    fitresult, meta_state, report = fit(tm, 0, X, y);
    history, _, state = meta_state;
    results2 = map(event -> last(event).measurement[1], history)
    @test results2 ≈ results
    @test fitresult.model == collect(r)[best_index]
    @test report.best_model == collect(r)[best_index]
    @test report.history == history
end

@testset_accelerated "accel. resampling" accel (exclude=[CPUThreads],) begin
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=[rms, l1], acceleration_resampling=accel)
    fitresult, meta_state, report = fit(tm, 0, X, y);
    history, _, state = meta_state;
    results3 = map(event -> last(event).measurement[1], history)
    @test results3 ≈ results
end

@testset_accelerated("under/over supply of models", accel,
                     (exclude=[CPUThreads],), begin
                     tm = TunedModel(model=first(r), tuning=Explicit(),
                                     range=r, measures=[rms, l1],
                                     acceleration=accel,
                                     resampling=CV(nfolds=2),
                                     n=4)
    mach = machine(tm, X, y)
    fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> last(event).measurement[1], history) ≈ results[1:4]

    tm.n=100
    @test_logs (:info, r"Only 12") fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> last(event).measurement[1], history) ≈ results
end)

end

true


text:
push!(LOAD_PATH,"../src/")

using Documenter
using QCDNUM
using Literate

# Generate examples notebooks
gen_content_dir = joinpath(@__DIR__, "src")

example_src = joinpath(@__DIR__, "..", "examples", "example.jl")
testsgns_src = joinpath(@__DIR__, "..", "examples", "testsgns.jl")
timing_src = joinpath(@__DIR__, "..", "examples", "timing.jl")
splint_src = joinpath(@__DIR__, "..", "examples", "splint.jl")

Literate.markdown(example_src, gen_content_dir, name="example")
Literate.markdown(testsgns_src, gen_content_dir, name="testsgns")
Literate.markdown(timing_src, gen_content_dir, name="timing")
Literate.markdown(splint_src, gen_content_dir, name="splint")

About = "Introduction" => "index.md"

Installation = "Installation" => "installation.md"

Examples = "Examples" => ["example.md", "testsgns.md",
                          "timing.md", "splint.md"]

Functions = "Available functions" => "functions.md"

Notebooks = "Notebook tutorial" => "notebook.md"

PAGES = [About, Installation, Notebooks, Examples, Functions]

makedocs(modules=[QCDNUM], sitename="QCDNUM.jl", pages=PAGES)

deploydocs(
    devbranch = "main",
    repo = "github.com/cescalara/QCDNUM.jl.git",
)


text:
module Day16

export get_inputs, get_solution1, get_solution2

using UnPack: @unpack
using InvertedIndices: Not

## Struct definitions
# Ticket rule
struct Rule
    field::String
    ranges::Vector{UnitRange}
end
function Rule(str)
    field, ranges = split(str, ':')
    ranges = parse_range.(split(ranges, " or "))
    return Rule(field, ranges)
end

## Input getting
function get_inputs()
    test_input1 = read_input("test_input1.txt")
    test_output1 = 71
    test_input2 = (read_input("test_input2.txt"), "")
    test_output2 = 12 * 11 * 13
    data = read_input("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_input(f_name)
    input = read(joinpath(@__DIR__, f_name), String)
    rules, my_ticket, nearby_tickets = split(input, r"\n\n|\r\n\r\n", keepempty=false)
    rules = Rule.(split(rules, r"\n|\r\n"))
    my_ticket = parse_ticket(split(my_ticket, r"\n|\r\n", limit=2)[2])
    nearby_tickets = parse_ticket.(split(nearby_tickets, r"\n|\r\n", keepempty=false)[2:end])
    return (; rules, my_ticket, nearby_tickets)
end

parse_range(str) = UnitRange(parse.(Int, split(str, '-'))...)

parse_ticket(str) = parse.(Int, split(str, ','))

## Solution functions
is_valid_under(number, rule) = any(number .∈ rule.ranges)
is_valid_under(ticket::AbstractArray, rules) = all(any(is_valid_under.(number, rules)) for number in ticket)
is_valid_under(rules) = x -> is_valid_under(x, rules)

drop_apply(f, x; dims) = dropdims(f(x; dims); dims)

# Part 1
function get_solution1(data)
    @unpack rules, nearby_tickets = data
    sum_invalids = 0
    for ticket in nearby_tickets, number in ticket
        if !any(is_valid_under.(number, rules))
            sum_invalids += number
        end
    end
    return sum_invalids
end

# Part 2
get_solution2(input::Tuple) = get_solution2(input...)

function get_solution2(data, word="departure")
    @unpack rules, my_ticket, nearby_tickets = data

    valid_tickets = filter(is_valid_under(rules), nearby_tickets)
    pushfirst!(valid_tickets, my_ticket)

    valid_matrix = [
        all(is_valid_under(ticket[i], rule) for ticket in valid_tickets)
        for rule in rules, i in eachindex(my_ticket)
    ]
    for i in sortperm(drop_apply(count, valid_matrix; dims=1))
        i_r = findfirst(valid_matrix[:,i])
        valid_matrix[i_r, Not(i)] .= false
    end
    rule_indices = findfirst.(eachcol(valid_matrix))

    return  prod(my_ticket[i] for (i, i_r) in enumerate(rule_indices) if occursin(word, rules[i_r].field))
end

end


text:
using SpecialExponentialFamilies
using Empirikos
using ApproxFun
using Random
using Intervals

expfamily = ExponentialFamily(basemeasure=Uniform(-4, 6), df=8, scale=true)

n = 10_000
α_smoothed_twin_tower =  [10; 0.0; 6.0; 9.0; -3.0; -2.0; -8.0; -1.0]

smoothed_twin_tower = expfamily(α_smoothed_twin_tower)

Random.seed!(1000)
Zs = rand(smoothed_twin_tower, n)
Zs_normal = StandardNormalSample.(Zs .+ randn(n))

_ef_fit = StatsBase.fit(SEF.PenalizedMLE(ef=expfamily, c0=1e-6), summarize(discretize(Zs_normal)))
_ef_fit.α_opt

ts = -3:0.1:3

lfsrs = Empirikos.PosteriorProbability.(StandardNormalSample.(ts), Intervals.Interval(0,nothing))
postmeans = Empirikos.PosteriorMean.(StandardNormalSample.(ts))

lfsr_cis = confint.(_ef_fit, lfsrs)
postmean_cis = confint.(_ef_fit, postmeans)

#using Plots
#pgfplotsx()
#plot(lfsrs, smoothed_twin_tower)
#plot!(ts, lfsr_cis)

#plot(postmeans, smoothed_twin_tower)
#plot!(ts, postmean_cis)


text:
using HttpCommon
using FactCheck
using HttpServer

facts("HttpServer utility functions:") do
    context("`write` correctly writes data response") do
        response = Response(200, "Hello World!")
        buf = IOBuffer();
        HttpServer.write(buf, response)
        response_string = String(take!(buf))
        vals = split(response_string, "\r\n")
        grep(a::Array, k::AbstractString) = filter(x -> ismatch(Regex(k), x), a)[1]
        @fact grep(vals, "HTTP") --> "HTTP/1.1 200 OK "
        @fact grep(vals, "Server") --> "Server: Julia/$VERSION"
        # default to text/html
        @fact grep(vals, "Content-Type") --> "Content-Type: text/html; charset=utf-8"
        # skip date
        @fact grep(vals, "Content-Language") --> "Content-Language: en"
        @fact grep(vals, "Hello") --> "Hello World!"
    end
end

import Requests: get, text, statuscode

facts("HttpServer runs") do
    context("using HTTP protocol on 0.0.0.0:8000") do
        http = HttpHandler() do req::Request, res::Response
            res = Response( ismatch(r"^/hello/",req.resource) ? string("Hello ", split(req.resource,'/')[3], "!") : 404 )
            setcookie!(res, "sessionkey", "abc", Dict("Path"=>"/test", "Secure"=>""))
        end
        server = Server(http)
        @async run(server, 8000)
        sleep(1.0)

        ret = Requests.get("http://localhost:8000/hello/travis")

        @fact text(ret) --> "Hello travis!"
        @fact statuscode(ret) --> 200
        @fact haskey(ret.cookies, "sessionkey") --> true

        let cookie = ret.cookies["sessionkey"]
            @fact cookie.value --> "abc"
            @fact cookie.attrs["Path"] --> "/test"
            @fact haskey(cookie.attrs, "Secure") --> true
        end

        ret = Requests.get("http://localhost:8000/bad")
        @fact text(ret) --> ""
        @fact statuscode(ret) --> 404
        close(server)
    end

    context("Rerun test using HTTP protocol on 0.0.0.0:8000 after closing") do
        http = HttpHandler() do req::Request, res::Response
            res = Response( ismatch(r"^/hello/",req.resource) ? string("Hello ", split(req.resource,'/')[3], "!") : 404 )
            setcookie!(res, "sessionkey", "abc", Dict("Path"=>"/test", "Secure"=>""))
        end
        server = Server(http)
        @async run(server, 8000)
        sleep(1.0)

        ret = Requests.get("http://localhost:8000/hello/travis")

        @fact text(ret) --> "Hello travis!"
        @fact statuscode(ret) --> 200
        @fact haskey(ret.cookies, "sessionkey") --> true

        let cookie = ret.cookies["sessionkey"]
            @fact cookie.value --> "abc"
            @fact cookie.attrs["Path"] --> "/test"
            @fact haskey(cookie.attrs, "Secure") --> true
        end

        ret = Requests.get("http://localhost:8000/bad")
        @fact text(ret) --> ""
        @fact statuscode(ret) --> 404
        close(server)
    end

    context("using HTTP protocol on 127.0.0.1:8001") do
        http = HttpHandler() do req::Request, res::Response
            Response( ismatch(r"^/hello/",req.resource) ? string("Hello ", split(req.resource,'/')[3], "!") : 404 )
        end
        server = Server(http)
        @async run(server, host=ip"127.0.0.1", port=8001)
        sleep(1.0)

        ret = Requests.get("http://127.0.0.1:8001/hello/travis")
        @fact text(ret) --> "Hello travis!"
        @fact statuscode(ret) --> 200
        close(server)
    end

    context("Rerun test using HTTP protocol on 127.0.0.1:8001 after closing") do
        http = HttpHandler() do req::Request, res::Response
            Response( ismatch(r"^/hello/",req.resource) ? string("Hello ", split(req.resource,'/')[3], "!") : 404 )
        end
        server = Server(http)
        @async run(server, host=ip"127.0.0.1", port=8001)
        sleep(1.0)

        ret = Requests.get("http://127.0.0.1:8001/hello/travis")
        @fact text(ret) --> "Hello travis!"
        @fact statuscode(ret) --> 200
        close(server)
    end

    context("Testing HTTPS on port 8002") do
        http = HttpHandler() do req, res
            Response("hello")
        end
        server = Server(http)
        cert = MbedTLS.crt_parse_file(joinpath(dirname(@__FILE__),"cert.pem"))
        key = MbedTLS.parse_keyfile(joinpath(dirname(@__FILE__),"key.pem"))
        @async run(server, port=8002, ssl=(cert, key))
        sleep(1.0)
        client_tls_conf = Requests.TLS_VERIFY
        MbedTLS.ca_chain!(client_tls_conf, cert)
        ret = Requests.get("https://localhost:8002", tls_conf=client_tls_conf)
        @fact text(ret) --> "hello"
        close(server)
    end

    context("Rerun test of HTTPS on port 8002 after closing") do
        http = HttpHandler() do req, res
            Response("hello")
        end
        server = Server(http)
        cert = MbedTLS.crt_parse_file(joinpath(dirname(@__FILE__),"cert.pem"))
        key = MbedTLS.parse_keyfile(joinpath(dirname(@__FILE__),"key.pem"))
        @async run(server, port=8002, ssl=(cert, key))
        sleep(1.0)
        client_tls_conf = Requests.TLS_VERIFY
        MbedTLS.ca_chain!(client_tls_conf, cert)
        ret = Requests.get("https://localhost:8002", tls_conf=client_tls_conf)
        @fact text(ret) --> "hello"
        close(server)
    end

    # Issue #111
    context("Parse HTTP headers") do
          http = HttpHandler() do req::Request, res::Response
              Response(req.headers["Content-Type"])
          end
          server = Server(http)
          @async run(server, 8000)
          sleep(1.0)

          ret = Requests.post("http://localhost:8000/", data = "√", headers = Dict("Content-Type" => "text/plain"))
          @fact text(ret) --> "text/plain"
          close(server)
    end
end


text: const ASSET_FINGERPRINT = "8d9151df5a4a5fafb268"

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
# example assuming * operation is always binary

# ENV["JULIA_DEBUG"] = Metatheory

using Metatheory

include("numberfold.jl")

comm_monoid = @theory begin
    ~a * ~b --> ~b * ~a
    ~a * 1 --> ~a
    ~a * (~b * ~c) --> (~a * ~b) * ~c
end

G = EGraph(:(3 * 4))
analyze!(G, NumberFold)

# exit(0)

@testset "Basic Constant Folding Example - Commutative Monoid" begin
    @test (true == @areequalg G comm_monoid 3 * 4 12)

    @test (true == @areequalg G comm_monoid 3 * 4 12 4 * 3  6 * 2)
end

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
    ex = :(a * 3 * b * 4)
    G = EGraph(ex)
    analyze!(G, NumberFold)
    addexpr!(G, :(12 * a))
    println(saturate!(G, comm_monoid))
    display(G.classes); println()
    @test (true == @areequalg G comm_monoid (12 * a) * b ((6 * 2) * b) * a)
    @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12 * a) * b ((6 * 2) * b) * a)
end

@testset "Basic Constant Folding Example - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    # addexpr!(G, 12)
    saturate!(G, comm_monoid)
    addexpr!(G, :(a * 2))
    analyze!(G, NumberFold)
    saturate!(G, comm_monoid)

    # display(G.classes); println()
    # println(G.root)
    # display(G.analyses[1].data); println()

    @test (true == areequal(G, comm_monoid, :(3 * 4), 12, :(4 * 3), :(6 * 2)))

    ex = :(a * 3 * b * 4)
    G = EGraph(ex)
    analyze!(G, NumberFold)
    params = SaturationParams(timeout=15)
    @test areequal(G, comm_monoid, :((3 * a) * (4 * b)), :((12 * a) * b),
        :(((6 * 2) * b) * a); params=params)
end

@testset "Infinite Loops analysis" begin
    boson = @theory begin
        1 * ~x --> ~x
    end

    G = EGraph(:(1 * x))
    params = SaturationParams(timeout=100)
    saturate!(G, boson, params)
    ex = extract!(G, astsize)

    # println(ex)

    using Metatheory.EGraphs
    boson = @theory begin
        (:c * :cdag) --> :cdag * :c + 1
        ~a * (~b + ~c) --> (~a * ~b) + (~a * ~c)
        (~b + ~c) * ~a --> (~b * ~a) + (~c * ~a)
        # 1 * x => x
        (~a * ~b) * ~c --> ~a * (~b * ~c)
        ~a * (~b * ~c) --> (~a * ~b) * ~c
    end

    g = EGraph(:(c * c * cdag * cdag))
    saturate!(g, boson)
    ex = extract!(g, astsize_inv)

    # println(ex)
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x'(t) = lm - d * x(t) - beta * x(t) * v(t),
  y'(t) = beta * x(t) * v(t) - a * y(t),
  v'(t) = k * y(t) - u * v(t),
  w'(t) = c * z(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
  z'(t) = c * q * y(t) * w(t) - h * z(t),
  y1(t) = w(t),
  y2(t) = z(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
include("structs.jl")
include("thermo.jl")

## Avoid typos when working with partial derivatives

## NW--N--NE
##     |
## W---C---E
##     |
## SW--S--SE

function cnswe(arr::AbstractMatrix{T},
        i::Int, j::Int) where T<:AbstractFloat
    c = arr[i, j]
    n = arr[i, j + 1]
    s = arr[i, j - 1]
    w = arr[i - 1, j]
    e = arr[i + 1, j]
    return c, n, s, w, e
end

function f_phase(s::AbstractFloat)::AbstractFloat
    return s^2
end

# ------------------ Continuity Equation ----------------- #

## Continuity equation to calculate ρᵢ
## φ ⋅∂ρᵢ / ∂t + div(ρᵢv⃗ᵢ) = 0, where ρᵢ = mᵢ / V
## (We can add a right side to this if we inject matter.)
function continuity_equation!(;
        ρnext::AbstractMatrix{T},
        ρwork::AbstractMatrix{T},
        F::AbstractMatrix{T},
        p::Parameters) where T<:AbstractFloat
    for j = 2 : p.ny - 1
        for i = 2 : p.nx - 1
            ρc, ρn, ρs, ρw, ρe = cnswe(ρwork, i, j)
            ρnext[i, j] = ρc + p.Δt  * F[i, j]
        end
    end
    @debug ρnext ρwork
end

function continuity_equation_fill_F!(;
        ρwork::AbstractMatrix{T},
        uwork::AbstractMatrix{T},
        vwork::AbstractMatrix{T},
        F::AbstractMatrix{T},
        p::Parameters) where T<:AbstractFloat
    for j = 2 : p.ny - 1
        for i = 2 : p.nx - 1
            ρc, ρn, ρs, ρw, ρe = cnswe(ρwork, i, j)
            uc, un, us, uw, ue = cnswe(uwork, i, j)
            vc, vn, vs, vw, ve = cnswe(vwork, i, j)

            F[i, j] = - 1 / p.φ *
                (uc * (ρe - ρw) / (2 * p.Δx) +
                 vc * (ρn - ρs) / (2 * p.Δy) +
                 ρc * (ue - uw) / (2 * p.Δx) +
                 ρc * (vn - vs) / (2 * p.Δy))
        end
    end
    @debug ρnext ρwork
end

# -------------------- Predictor ------------------------- #
function continuity_equation!(syswork::System,
        sysnext::System, p::Parameters)
    for k = 1 : 2
        @debug "Component $k continuity eqn"
        @views continuity_equation_fill_F!(
                                    ρwork=syswork.ρ[:, :, k],
                                    uwork=syswork.u[:, :, k],
                                    vwork=syswork.v[:, :, k],
                                    F=syswork.F[:, :, k], p=p)
        @views continuity_equation!(ρnext=sysnext.ρ[:, :, k],
                                    ρwork=syswork.ρ[:, :, k],
                                    F=syswork.F[:, :, k], p=p)
    end
    @debug "Saturation " sysnext.s
end
# -------------------------------------------------------- #

# -------------------- Corrector ------------------------- #
function continuity_equation!(syswork::System,
        sysnext_pred::System, sysnext::System,
        p::Parameters)
    for k = 1 : 2
        @debug "Component $k continuity eqn"
        @views continuity_equation_fill_F!(
                                    ρwork=sysnext_pred.ρ[:, :, k],
                                    uwork=sysnext_pred.u[:, :, k],
                                    vwork=sysnext_pred.v[:, :, k],
                                    F=sysnext_pred.F[:, :, k], p=p)
        @views continuity_equation!(ρnext=sysnext.ρ[:, :, k],
                                    ρwork=syswork.ρ[:, :, k],
                                    F=(syswork.F[:, :, k] .+
                                       sysnext_pred.F[:, :, k]) ./2,
                                    p=p)
    end
    @debug "Saturation " sysnext.s
end
# -------------------------------------------------------- #

# ---------------------- Darcy --------------------------- #

## Darcy's law
## v⃗ = -μ⁻¹ K̂⋅f(s) ⋅∇P

function darcy!(;
        u::AbstractMatrix{T}, v::AbstractMatrix{T},
        P::AbstractMatrix{T}, s::AbstractMatrix{T},
        p::Parameters, k::Integer) where T<:AbstractFloat
    for j = 2 : p.ny - 1
        for i = 2 : p.nx - 1
            ## можно через (ρwork + ρ) / 2 -> 2-го порядка
            Pc, Pn, Ps, Pw, Pe = cnswe(P, i, j)

            u[i, j] = -p.K / p.μ[k] * f_phase(s[i, j]) *
                (Pe - Pw) / (2 * p.Δx)
                v[i, j] = -p.K / p.μ[k] * f_phase(s[i, j]) *
                (Pn - Ps) / (2 * p.Δy)
        end
    end
end

function darcy!(sys::System, p::Parameters)
    ## Same problem with type-stability as boundary_velocity
    ## https://discourse.julialang.org/t/how-to-make-type-stable-passing-an-element-of-a-vector-of-functions/60068
    # fgas(s) = s^2
    # fliq(s) = (1 .- s)^2
    # f = (fgas, fliq)

    # f(s::AbstractFloat)::AbstractFloat = s^2

    for k = 1 : 2
      @views darcy!(u=sys.u[:, :, k],
                    v=sys.v[:, :, k], P=sys.P,
                    s=sys.s[:, :, k], p=p, k=k)
    end

    ## This is type stable
    ## (but I'm still getting runtime dispatch in profiler)
    # fgas(s) = s^2
    # fliq(s) = (1 .- s)^2

    # @views darcy!(f=fgas, u=sys.u[:, :, 1],
    #               v=sys.v[:, :, 1], P=sys.P,
    #               s=sys.s, p=p, k=1)

    # @views darcy!(f=fliq, u=sys.u[:, :, 2],
    #               v=sys.v[:, :, 2], P=sys.P,
    #               s=sys.s, p=p, k=2)
end
# -------------------------------------------------------- #

# ---------------------- BC & IC ------------------------- #

function initial_saturation!(
        ρ̂₁::AbstractMatrix{T},
        ρ̂₂::AbstractMatrix{T},
        s::AbstractArray{T, 3}, p) where T<:AbstractFloat

    @. @views s[:, :, 1] = (ρ̂₁[:, :] /
                         ρ̂₂[:, :] * p.M[2] * (1 - p.ψ₀) /
                         p.M[1] / p.ψ₀ + 1)^(-1)

    @. @views s[:, :, 2] = 1 - s[:, :, 1]
end

function initial_saturation!(sys::System, p::Parameters)
    ρ̂₁ = density.(ideal_gas, sys.P)
    ρ̂₂ = density.(tait_C₅H₁₂, sys.P)
    initial_saturation!(ρ̂₁, ρ̂₂, sys.s, p)
end

function apply_bc!(; grid_ghost::AbstractMatrix{T},
        grid_inside::AbstractMatrix{T}, a::AbstractMatrix{T},
        b::AbstractMatrix{T}, rhs::AbstractMatrix{T},
        p::Parameters, Δh::T) where T<:AbstractFloat
    @views @. grid_ghost = (2 * Δh * rhs + grid_inside *
                                   (2 * b - a * Δh)) /
                                  (a * Δh + 2 * b)
end

function apply_bc!(grid::Array{T, 3}, bc::BoundaryCondition{T},
        rhs::Array{T, 3}, p::Parameters) where T<:AbstractFloat
    ## left
    @views apply_bc!(grid_ghost  = grid[1, :, :],
              grid_inside = grid[2, :, :],
              a   = bc.a[:, 1, :],
              b   = bc.b[:, 1, :],
              rhs =  rhs[:, 1, :], p=p, Δh = p.Δx)

    ## right
    @views apply_bc!(grid_ghost  = grid[nx,     :, :],
              grid_inside = grid[nx - 1, :, :],
              a   = bc.a[:, 2, :],
              b   = bc.b[:, 2, :],
              rhs =  rhs[:, 2, :], p=p, Δh = p.Δx)

    ## top
    @views apply_bc!(grid_ghost  = grid[:, ny,     :],
              grid_inside = grid[:, ny - 1, :],
              a   = bc.a[:, 3, :],
              b   = bc.b[:, 3, :],
              rhs =  rhs[:, 3, :], p=p, Δh = p.Δy)

    ## bottom
    @views apply_bc!(grid_ghost  = grid[:, 1, :],
              grid_inside = grid[:, 2, :],
              a   = bc.a[:, 4, :],
              b   = bc.b[:, 4, :],
              rhs =  rhs[:, 4, :], p=p, Δh = p.Δy)
end

function update_RHS_saturation!(; s_RHS::AbstractMatrix{T},
        P::Array{T, 1}, p::Parameters,
        t::Integer) where T<:AbstractFloat
    @views ρ̂₁_bound    = density.(ideal_gas,  P)
    @views ρ̂₂_bound    = density.(tait_C₅H₁₂, P)

    if t < p.relax_steps
        ψ = (p.ψ - p.ψ₀) * t / relax_steps + p.ψ₀
    else
        ψ = p.ψ
    end

    @views s_RHS[:, 1] = saturation.(ρ̂₁_bound, ρ̂₂_bound, ψ,
                                         p.M[1], p.M[2])
    @views @. s_RHS[:, 2] = 1 - s_RHS[:, 1]
end

function update_RHS_saturation!(s_RHS::Array{T, 3}, P::Array{T, 2},
        p::Parameters, t::Integer) where T<:AbstractFloat
    ## bottom inlet
    @views update_RHS_saturation!(s_RHS = s_RHS[p.inlet, 4, :],
                           P = fill(p.Pin, last(collect(p.inlet))),
                           p = p, t = t)
    ## everywhere else we don't have to do anything because
    ## ∂s / ∂n = 0
end

function get_bound(h::AbstractArray{T, 2},
        p::Parameters) where T<:AbstractFloat
    @views h_left  = (h[1, :]    + h[2, :])        / 2
    @views h_right = (h[p.nx, :] + h[p.nx - 1, :]) / 2
    @views h_top   = (h[:, p.ny] + h[:, p.ny - 1]) / 2
    @views h_bott  = (h[:, 1]    + h[:, 2])        / 2
    return hcat(h_left, h_right, h_top, h_bott)
end

function update_RHS_density!(ρ_RHS::Array{T, 3},
        s::Array{T, 3}, P::Array{T, 2},
        p::Parameters) where T<:AbstractFloat
    @views P_bound  = get_bound(P, p)
    @views s₁_bound = get_bound(s[:, :, 1], p)
    @views s₂_bound = 1 .- s₁_bound

    @views ρ̂₁_bound = density.(ideal_gas,  P_bound)
    @views ρ̂₂_bound = density.(tait_C₅H₁₂, P_bound)
    @views @. ρ_RHS[:, :, 1] = s₁_bound * ρ̂₁_bound
    @views @. ρ_RHS[:, :, 2] = s₂_bound * ρ̂₂_bound
end

function update_RHS_pressure!(P_RHS::Array{T, 3},
        p::Parameters, t::Integer) where T<:AbstractFloat
    if t < p.relax_steps
        P_inlet = (p.Pin - p.P₀) * t / p.relax_steps + p.P₀
    else
        P_inlet = p.Pin
    end

    ## bottom inlet
    @views @. P_RHS[p.inlet, 4, 1] = P_inlet

    ## top outlet
    @views @. P_RHS[:, 3, 1]       = p.Pout

    ## at walls ∂P / ∂n = 0 so we don't have to change anything,
    ## since P_RHS = 0
end

function update_RHS_velocity!(v_RHS::Array{T, 3},
        P::Array{T, 2}, s::Array{T, 3},
        p::Parameters) where T<:AbstractFloat
    @views s₁_bound = get_bound(s[:, :, 1], p)
    @views s₂_bound = 1 .- s₁_bound
    @views s_bound = reshape(hcat(s₁_bound, s₂_bound),
                      (max(p.nx, p.ny), 4, 2))

    ## bottom inlet
    for k = 1 : 2
        @. @views v_RHS[p.inlet, 4, k] =  -p.K / p.μ[k] *
                                  f_phase(s_bound[p.inlet, 4, k]) *
                                  (P[p.inlet, 2] -
                                   P[p.inlet, 1]) / p.Δy

        ## top outlet
        @. @views v_RHS[:, 3, k]       =  -p.K / p.μ[k] *
                                    f_phase(s_bound[:, 3, k]) *
                                    (P[:, ny] -
                                     P[:, ny - 1]) / p.Δy
    end
end

function apply_bc_pres_sat_den!(sys::System, p::Parameters,
        t::Integer)
    ## ∂P / ∂n = 0 at walls
    ## P = Pin     at inlet
    ## P = Pout    at outlet
    update_RHS_pressure!(sys.bc_RHS.P, p, t)
    P_tmp = zeros(p.nx, p.ny, 2)
    P_tmp[:, :, 1] = sys.P
    apply_bc!(P_tmp, p.bc.P, sys.bc_RHS.P, p)
    sys.P = P_tmp[:, :, 1]

    ## s_inlet = s_inlet(ρ̂₁(T, P), ρ̂₂(P, T), P, T)
    ## ∂s / ∂n = 0 everywhere else
    update_RHS_saturation!(sys.bc_RHS.s, sys.P, p, t)
    apply_bc!(sys.s, p.bc.s, sys.bc_RHS.s, p)

    ## ρ_k_bound = ρ̂_k_bound(T, P) * s_k_bound(T, P)
    update_RHS_density!(sys.bc_RHS.ρ, sys.s, sys.P, p)
    apply_bc!(sys.ρ, p.bc.ρ, sys.bc_RHS.ρ, p)
end

function apply_bc_velocity!(sys::System, p::Parameters)
    ## Darcy at inlet and outlet
    ## v_n at walls
    update_RHS_velocity!(sys.bc_RHS.v, sys.P, sys.s, p)
    apply_bc!(sys.u, p.bc.u, sys.bc_RHS.u, p)
    apply_bc!(sys.v, p.bc.v, sys.bc_RHS.v, p)
end

function initial_conditions!(sys::System, p::Parameters)
    ## ∂P / ∂n = 0 at walls
    ## P = Pin     at inlet
    ## P = Pout    at outlet
    update_RHS_pressure!(sys.bc_RHS.P, p, 0)
    P_tmp = zeros(p.nx, p.ny, 2)
    P_tmp[:, :, 1] = sys.P
    apply_bc!(P_tmp, p.bc.P, sys.bc_RHS.P, p)
    sys.P = P_tmp[:, :, 1]

    initial_saturation!(sys, p)

    ## s_inlet = s_inlet(ρ̂₁(T, P), ρ̂₂(P, T), P, T)
    update_RHS_saturation!(sys.bc_RHS.s, sys.P, p, 0)
    apply_bc!(sys.s, p.bc.s, sys.bc_RHS.s, p)

    density!(sys)
    update_RHS_density!(sys.bc_RHS.ρ, sys.s, sys.P, p)
    apply_bc!(sys.ρ, p.bc.ρ, sys.bc_RHS.ρ, p)

    ## Darcy at inlet and outlet
    update_RHS_velocity!(sys.bc_RHS.v, sys.P, sys.s, p)
    apply_bc!(sys.u, p.bc.u, sys.bc_RHS.u, p)
    apply_bc!(sys.v, p.bc.v, sys.bc_RHS.v, p)
end

function initial_RHS(p::Parameters)
    u_RHS = zeros(max(p.nx, p.ny), 4, 2)
    v_RHS = zeros(max(p.nx, p.ny), 4, 2)
    ρ_RHS = zeros(max(p.nx, p.ny), 4, 2)
    s_RHS = zeros(max(p.nx, p.ny), 4, 2)
    P_RHS = zeros(max(p.nx, p.ny), 4, 2)
    return BoundaryRHS(u_RHS, v_RHS, ρ_RHS, s_RHS, P_RHS)
end
# -------------------------------------------------------- #

# -------------------- Flux ------------------------------ #
function get_boundary(a::Array{T, 2},
        j::Integer) where T<:AbstractFloat
    @views a_bound = (a[:, j] + a[:, j + 1]) / 2
    return a_bound
end

function get_flux(ρ::Array{T, 2}, s::Array{T, 2}, v::Array{T, 2},
        j::Integer) where T<:AbstractFloat
    ρ_bound = get_boundary(ρ, j)
    s_bound = get_boundary(s, j)
    v_bound = get_boundary(v, j)
    return ρ_bound ./ s_bound .* v_bound
end

function get_total_flux(sys::System, p::Parameters)
    flux_init  = zeros(2, p.nx)
    flux_final = zeros(2, p.nx)
    for k = 1 : 2
        flux_init[k, :]  = get_flux(sys.ρ[:, :, k],
                                    sys.s[:, :, k],
                                    sys.v[:, :, k], 1)
        flux_final[k, :] = get_flux(sys.ρ[:, :, k],
                                    sys.s[:, :, k],
                                    sys.v[:, :, k], p.ny - 1)
    end
    flux_init_total  = sum(flux_init,  dims = 2)
    flux_final_total = sum(flux_final, dims = 2)
    return flux_init_total, flux_final_total
end

function get_content_inside(sys::System, p::Parameters)
    content = zeros(2)
    @views content[1] = sum(sys.ρ[:, :, 1] ./ sys.s[:, :, 1])
    @views content[2] = sum(sys.ρ[:, :, 2] ./ sys.s[:, :, 2])
    return content
end
# -------------------------------------------------------- #

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
__precompile__()

module WriteVTK

# All the code is based on the VTK file specification [1], plus some
# undocumented stuff found around the internet...
# [1] http://www.vtk.org/VTK/img/file-formats.pdf

export VTKCellTypes, VTKCellType
export MeshCell
export vtk_grid, vtk_save, vtk_point_data, vtk_cell_data
export vtk_multiblock
export paraview_collection, collection_add_timestep

using LightXML
using BufferedStreams: BufferedOutputStream, EmptyStreamSource
using Libz: ZlibDeflateOutputStream
import Compat.UTF8String

import Base: close, isopen

# Cell type definitions as in vtkCellType.h
include("VTKCellTypes.jl")

if VERSION >= v"0.5"
    Base.@deprecate_binding VTKCellType VTKCellTypes
else
    const VTKCellType = VTKCellTypes
end

## Constants ##
const COMPRESSION_LEVEL = 6
const IS_LITTLE_ENDIAN = (ENDIAN_BOM == 0x04030201)  # see the documentation for ENDIAN_BOM

## Types ##
abstract VTKFile

typealias DataBuffer BufferedOutputStream{EmptyStreamSource}

immutable DatasetFile <: VTKFile
    xdoc::XMLDocument
    path::UTF8String
    gridType_str::UTF8String
    Npts::Int           # Number of grid points.
    Ncls::Int           # Number of cells.
    compressed::Bool    # Data is compressed?
    appended::Bool      # Data is appended? (otherwise it's written inline, base64-encoded)
    buf::DataBuffer     # Buffer with appended data.
    function DatasetFile(xdoc, path, gridType_str, Npts, Ncls,
                         compressed, appended)
        buf = BufferedOutputStream() :: DataBuffer
        if !appended  # in this case we don't need a buffer
            close(buf)
        end
        return new(xdoc, path, gridType_str, Npts, Ncls,
                   compressed, appended, buf)
    end
end

immutable MultiblockFile <: VTKFile
    xdoc::XMLDocument
    path::UTF8String
    blocks::Vector{VTKFile}
    # Constructor.
    MultiblockFile(xdoc, path) = new(xdoc, path, VTKFile[])
end

immutable CollectionFile <: VTKFile
    xdoc::XMLDocument
    path::UTF8String
    timeSteps::Vector{VTKFile}
    # Constructor.
    CollectionFile(xdoc, path) = new(xdoc, path, VTKFile[])
end

# Cells in unstructured meshes.
immutable MeshCell
    ctype::VTKCellTypes.VTKCellType           # cell type identifier (see VTKCellTypes.jl)
    connectivity::Vector{Int32}  # indices of points (one-based, following the convention in Julia)

    function MeshCell{T<:Integer}(ctype::VTKCellTypes.VTKCellType,connectivity::Vector{T})
        if length(connectivity) == ctype.nodes || ctype.nodes == -1
            new(ctype,connectivity)
        else
            throw(ArgumentError("Wrong number of nodes in connectivity vector."))
        end
    end
end

close(vtk::VTKFile) = free(vtk.xdoc)
isopen(vtk::VTKFile) = (vtk.xdoc.ptr != C_NULL)

# Multiblock-specific functions and types.
include("gridtypes/multiblock.jl")
include("gridtypes/ParaviewCollection.jl")

# Grid-specific functions and types.
include("gridtypes/structured.jl")
include("gridtypes/unstructured.jl")
include("gridtypes/rectilinear.jl")
include("gridtypes/imagedata.jl")

# Common functions.
include("gridtypes/common.jl")

end


text:
# configurations
# --------------

# TODO more configurations, e.g. ignore user-specified modules and such
"""
Configurations for JET analysis.
These configurations will be active for all the entries.

---
- `strict_condition_check::Bool = false` \\
  Enables strict condition check.
  JET reports an error if a condition expression type is "non-boolean". In a case when
    the condition type is `Union`, JET will report if either of union split case is
    non-boolean type, but this can lead to lots of false positive error reports when the
    code is not well-typed, because Julia `Base` defines generic functions that are commonly
    used at a conditional context but also may return "non-boolean" values, e.g.:
    - `!(::Function) -> Function`
    - `!(::Missing) -> Missing`
    - `==(::Missing, ::Any) -> Missing`
    - `==(::Any, ::Missing) -> Missing`
    and thus loosely-typed conditional expression often becomes e.g. `Union{Bool, Missing}`,
    and consequently JET will report it as "non-boolean" type
    (NOTE: in Julia `Missing` is certainly not valid conditional type).
  If this configuration is set to `false`, JET enables an heuristic to avoid those false
    positive error reports and won't report an error if a condition expression type is
    `Union` and either of its union split case is `Function` or `Missing`.

  The effect of this configuration can be described with the following examples:

  * with `strict_condition_check::Bool = false` (default)
    ```julia
    julia> test_f() = Dict('a' => 1, :b => 2) # ::Dict{Any,Int}
    test_f (generic function with 1 method)

    julia> @report_call test_f()
    No errors !
    Dict{Any, Int64}
    ```

  * with `strict_condition_check::Bool = true`
    ```julia
    julia> test_f() = Dict('a' => 1, :b => 2) # ::Dict{Any,Int}
    test_f (generic function with 1 method)

    julia> @report_call strict_condition_check = true test_f()
    ═════ 1 possible error found ═════
    ┌ @ REPL[2]:1 Main.Dict(Main.=>('a', 1), Main.=>(:b, 2))
    │┌ @ dict.jl:125 Base.Dict(ps)
    ││┌ @ dict.jl:129 Base.dict_with_eltype(#308, kv, Base.eltype(kv))
    │││┌ @ abstractdict.jl:539 Base.grow_to!(Base.dict_with_eltype(DT_apply, _5), kv)
    ││││┌ @ dict.jl:145 Base.grow_to!(dest2, itr, st)
    │││││┌ @ dict.jl:159 Base.setindex!(new, v, k)
    ││││││┌ @ dict.jl:383 Base.ht_keyindex2!(h, key)
    │││││││┌ @ dict.jl:328 goto %35 if not Base.isequal(key, Base.getindex(keys, index))
    ││││││││ for 1 of union split cases, non-boolean (Missing) used in boolean context: goto %35 if not Base.isequal(key::Symbol, Base.getindex(keys::Vector{Any}, index::Int64)::Any)::Union{Missing, Bool}
    │││││││└───────────────
    Dict{Any, Int64}
    ```
---
- `ignore_native_remarks::Bool = true` \\
  If `true`, JET won't construct nor cache reports of "native remarks", which may speed up analysis time.
  "Native remarks" are information that Julia's native compiler emits about how type inference routine goes,
  and those remarks are less interesting in term of "error checking", so JET ignores them by default.
"""
struct JETAnalysisParams
    strict_condition_check::Bool
    ignore_native_remarks::Bool
    @jetconfigurable JETAnalysisParams(; strict_condition_check::Bool = false,
                                         ignore_native_remarks::Bool  = true,
                                         ) =
        return new(strict_condition_check,
                   ignore_native_remarks,
                   )
end

"""
Configurations for Julia's native type inference routine.
These configurations will be active for all the entries.

You can specify all the keyword parameters of `Core.Compiler.InferenceParams`, e.g.
  `max_methods::Int = 3`, `union_splitting::Int = 4`.
Listed here are selections of those parameters that can have a potent influence on JET analysis.

---
- `ipo_constant_propagation::Bool = true` \\
  Enables constant propagation in abstract interpretation.
  It is _**highly**_ recommended that you keep this configuration `true` to get reasonable analysis,
  because constant propagation can cut off lots of false positive errorenous code paths and
  thus lead to more accurate and useful analysis results.
---
- `aggressive_constant_propagation::Bool = true` \\
  If `true`, JET will try to do constant propagation more "aggressively".
  As explained above, it can lead to more accurate analysis, but also lead to worse analysis
  performance at the cost of that.
---
- `unoptimize_throw_blocks::Bool = false` \\
  Turn this on to skip analysis on code blocks that will eventually lead to a `throw` call.
  This configuration may improve the analysis performance, but it's better to be turned off
    for JET analysis, because there may be other errors even in those code blocks.
"""
@jetconfigurable JETInferenceParams(; ipo_constant_propagation::Bool        = true,
                                      aggressive_constant_propagation::Bool = true,
                                      unoptimize_throw_blocks::Bool         = false,
                                      max_methods::Int                      = 3,
                                      union_splitting::Int                  = 4,
                                      apply_union_enum::Int                 = 8,
                                      tupletype_depth::Int                  = 3,
                                      tuple_splat::Int                      = 32,
                                      ) =
    return InferenceParams(; ipo_constant_propagation,
                             aggressive_constant_propagation,
                             unoptimize_throw_blocks,
                             max_methods,
                             union_splitting,
                             apply_union_enum,
                             tupletype_depth,
                             tuple_splat,
                             )

# here we just try to sync `OptimizationParams` with  `InferenceParams`
# (the same JET configurations will be passed on here)
JETOptimizationParams(; # inlining::Bool                = inlining_enabled(),
                        # inline_cost_threshold::Int    = 100,
                        # inline_nonleaf_penalty::Int   = 1000,
                        # inline_tupleret_bonus::Int    = 250,
                        # inline_error_path_cost::Int   = 20,
                        max_methods::Int              = 3,
                        tuple_splat::Int              = 32,
                        union_splitting::Int          = 4,
                        unoptimize_throw_blocks::Bool = false,
                        _jetconfigs...) =
    # NOTE we always disable inlining, because our current strategy to find undefined
    # local variable assumes un-inlined frames
    # TODO enable inlining to get better JET analysis performance ?
    # XXX but the self-profiling with `inlining = true` showed performance regression ...
    return OptimizationParams(; inlining = false,
                                # inline_cost_threshold,
                                # inline_nonleaf_penalty,
                                # inline_tupleret_bonus,
                                # inline_error_path_cost,
                                max_methods,
                                tuple_splat,
                                union_splitting,
                                unoptimize_throw_blocks,
                                )

const LOGGER_LEVEL_KEY     = :JET_LOGGER_LEVEL
const INFO_LOGGER_LEVEL    = 0
const DEBUG_LOGGER_LEVEL   = 1
const LOGGER_LEVELS        = Dict(INFO_LOGGER_LEVEL  => :info,
                                  DEBUG_LOGGER_LEVEL => :debug,
                                  )
const DEFAULT_LOGGER_LEVEL = INFO_LOGGER_LEVEL
const LOGGER_LEVELS_DESC   = let
    map(collect(LOGGER_LEVELS)) do (level, desc)
        if level == DEFAULT_LOGGER_LEVEL
            "`$level` (\"$desc\" level, default)"
        else
            "`$level` (\"$desc\" level)"
        end
    end |> Fix2(join, ", ")
end
get_logger_level(io::IO)   = get(io, LOGGER_LEVEL_KEY, DEFAULT_LOGGER_LEVEL)::Int

"""
Logging configurations for JET analysis.

---
- `toplevel_logger::Union{Nothing,IO} = nothing` \\
  If `IO` object is given, it will track JET's toplevel analysis.
  Logging level can be specified with `$(repr(LOGGER_LEVEL_KEY))` `IO` property.
  Currently supported logging levels are either of $(LOGGER_LEVELS_DESC).

  Examples:
  * logs into `stdout`
    ```julia
    julia> report_file(filename; toplevel_logger = stdout)
    ```
  * logs into `io::IOBuffer` with "debug" logger level
    ```julia
    julia> report_file(filename; toplevel_logger = IOContext(io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL));
    ```
---
- `inference_logger::Union{Nothing,IO} = nothing` \\
  If `IO` object is given, it will track JET's abstract interpretation routine.
  Logging level can be specified with `$(repr(LOGGER_LEVEL_KEY))` `IO` property.
  Currently supported logging levels are either of $(LOGGER_LEVELS_DESC).

  Examples:
  * logs into `stdout`
    ```julia
    analyze_call(f, args...; inference_logger = stdout)
    ```
  * logs into `io::IOBuffer` with "debug" logger level
    ```julia
    julia> analyze_call(f, args...; inference_logger = IOContext(io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL))
    ```
---
!!! tip
    Of course you can specify both `toplevel_logger` and `inference_logger` at the same time like below:
    ```julia
    julia> report_and_watch_file(filename;
                                 toplevel_logger = IOContext(logger_io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL),
                                 inference_logger = inference_io)
    ```
"""
struct JETLogger
    # logger to track toplevel interpretation
    toplevel_logger::Union{Nothing,IO}
    # logger to track abstract interpretation routine
    inference_logger::Union{Nothing,IO}
    @jetconfigurable function JETLogger(; toplevel_logger::Union{Nothing,IO} = nothing,
                                          inference_logger::Union{Nothing,IO} = nothing,
                                          )
        if isa(toplevel_logger, IO)
            @assert get_logger_level(toplevel_logger) in keys(LOGGER_LEVELS) "toplevel_logger's $LOGGER_LEVEL_KEY should be either of $LOGGER_LEVELS_DESC"
        end
        if isa(inference_logger, IO)
            @assert get_logger_level(inference_logger) in keys(LOGGER_LEVELS) "inference_logger's $LOGGER_LEVEL_KEY should be either of $LOGGER_LEVELS_DESC"
        end
        new(toplevel_logger, inference_logger)
    end
end

# XXX we need to consider world range ?
struct AnalysisResult
    linfo::MethodInstance
    argtypes::Vector{Any}
    overridden_by_const::CC.BitVector
    cache::Vector{InferenceErrorReportCache}
end

mutable struct JETInterpreter <: AbstractInterpreter
    ### native ###

    native::NativeInterpreter

    ### JET.jl specific ###

    ## general ##

    # key for the JET's global cache
    cache_key::UInt

    # reports found so far
    reports::Vector{InferenceErrorReport}

    # stashes `UncaughtExceptionReport`s that are not caught so far
    uncaught_exceptions::Vector{UncaughtExceptionReport}

    # keeps reports that should be updated when returning back the parent frame (i.e. the next time we get back to inter-procedural context)
    to_be_updated::Set{InferenceErrorReport}

    # keeps track of the current inference frame (needed for report cache reconstruction)
    current_frame::Union{Nothing,InferenceState}

    # local report cache for constant analysis
    cache::Vector{AnalysisResult}

    # configurations for JET analysis
    analysis_params::JETAnalysisParams

    # tracks there has been any error reported in `abstract_call_method`
    # that information will help JET decide it's worth to do constant prop' even if
    # `NativeInterpreter` doesn't find it useful (i.e. the return type can't be improved anymore)
    # NOTE this field is supposed to be computed at each call of `abstract_call_method`, and
    # then immediately reset within the next `abstract_call_method_with_const_args` call
    anyerror::Bool

    ## virtual toplevel execution ##

    # will be used in toplevel analysis (skip inference on actually interpreted statements)
    concretized::BitVector

    # virtual toplevel module
    toplevelmod::Module

    # toplevel modules concretized by JET (only active within sequential toplevel analysis)
    toplevelmods::Set{Module}

    # slots to represent toplevel global variables
    global_slots::Dict{Int,Symbol}

    ## debug ##

    logger::JETLogger

    # records depth of call stack
    depth::Int

    function JETInterpreter(native::NativeInterpreter, args...)
        @assert !native.opt_params.inlining "inlining should be disabled for JETInterpreter analysis"
        new(native, args...)
    end
end

# constructor for fresh analysis
@jetconfigurable function JETInterpreter(world::UInt     = get_world_counter();
                                         current_frame   = nothing,
                                         cache           = AnalysisResult[],
                                         analysis_params = nothing,
                                         inf_params      = nothing,
                                         opt_params      = nothing,
                                         concretized     = _CONCRETIZED,
                                         toplevelmod     = _TOPLEVELMOD,
                                         toplevelmods    = _TOPLEVELMODS,
                                         global_slots    = _GLOBAL_SLOTS,
                                         logger          = nothing,
                                         depth           = 0,
                                         jetconfigs...)
    isnothing(analysis_params) && (analysis_params = JETAnalysisParams(; jetconfigs...))
    isnothing(inf_params)      && (inf_params = JETInferenceParams(; jetconfigs...))
    isnothing(opt_params)      && (opt_params = JETOptimizationParams(; jetconfigs...))
    isnothing(logger)          && (logger = JETLogger(; jetconfigs...))

    cache_key = gen_cache_key(analysis_params, inf_params)
    haskey(JET_REPORT_CACHE, cache_key) || (JET_REPORT_CACHE[cache_key] = IdDict())
    haskey(JET_CODE_CACHE, cache_key) || (JET_CODE_CACHE[cache_key] = IdDict())

    return JETInterpreter(NativeInterpreter(world; inf_params, opt_params),
                          cache_key,
                          InferenceErrorReport[],
                          UncaughtExceptionReport[],
                          Set{InferenceErrorReport}(),
                          current_frame,
                          cache,
                          analysis_params,
                          false,
                          concretized,
                          toplevelmod,
                          toplevelmods,
                          global_slots,
                          logger,
                          depth,
                          )
end

# dummies for non-toplevel analysis
module __toplevelmod__ end
const _CONCRETIZED  = BitVector()
const _TOPLEVELMOD  = __toplevelmod__
const _TOPLEVELMODS = Set{Module}()
const _GLOBAL_SLOTS = Dict{Int,Symbol}()

# constructor for sequential toplevel JET analysis
function JETInterpreter(interp::JETInterpreter, concretized, toplevelmod)
    return JETInterpreter(# update world age to take in newly added methods defined
                          # in a previous toplevel interpretation performed by `ConcreteInterpreter`
                          get_world_counter();
                          analysis_params = JETAnalysisParams(interp),
                          inf_params      = InferenceParams(interp),
                          opt_params      = OptimizationParams(interp),
                          concretized     = concretized, # or construct partial `CodeInfo` from remaining abstract statements ?
                          toplevelmod     = toplevelmod,
                          toplevelmods    = push!(interp.toplevelmods, toplevelmod),
                          logger          = JETLogger(interp),
                          )
end

# constructor to do additional JET analysis in the middle of parent (non-toplevel) interpretation
function JETInterpreter(interp::JETInterpreter)
    return JETInterpreter(get_world_counter(interp);
                          current_frame   = interp.current_frame,
                          cache           = interp.cache,
                          analysis_params = JETAnalysisParams(interp),
                          inf_params      = InferenceParams(interp),
                          opt_params      = OptimizationParams(interp),
                          logger          = JETLogger(interp),
                          depth           = interp.depth,
                          )
end

function Base.show(io::IO, interp::JETInterpreter)
    rn = length(interp.reports)
    en = length(interp.uncaught_exceptions)
    print(io, "JETInterpreter with $(rn) reports and $(en) uncaught exceptions")
    frame = interp.current_frame
    if !isnothing(frame)
        print(io, " at ")
        show(io, frame.linfo)
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", interp::JETInterpreter) =
    return interp
Base.show(io::IO, ::MIME"application/prs.juno.inline", interp::NativeInterpreter) =
    return interp

# XXX: this should be upstreamed
function Base.show(io::IO, frame::InferenceState)
    print(io, "InfernceState for ")
    show(io, frame.linfo)
    print(io, " with currpc ", frame.currpc, '/', length(frame.src.code))
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", frame::InferenceState) =
    return frame

cache_key(interp::JETInterpreter) = interp.cache_key
function gen_cache_key(analysis_params::JETAnalysisParams, inf_params::InferenceParams)
    h = @static UInt === UInt64 ? 0xa49bd446c0a5d90e : 0xe45361ac
    h = hash(analysis_params, h)
    h = hash(inf_params, h)
    return h
end

# AbstractInterpreter API
# -----------------------

CC.InferenceParams(interp::JETInterpreter)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::JETInterpreter) = OptimizationParams(interp.native)
CC.get_world_counter(interp::JETInterpreter)  = get_world_counter(interp.native)

# JET only works for runtime inference
CC.lock_mi_inference(::JETInterpreter, ::MethodInstance) = nothing
CC.unlock_mi_inference(::JETInterpreter, ::MethodInstance) = nothing

# function CC.add_remark!(interp::JETInterpreter, sv::InferenceState, s::String)
#     JETAnalysisParams(interp).ignore_native_remarks && return
#     push!(interp.native_remarks, NativeRemark(interp, sv, s))
#     return
# end
CC.add_remark!(interp::JETInterpreter, sv::InferenceState, s::String) = return

CC.may_optimize(interp::JETInterpreter)      = true
CC.may_compress(interp::JETInterpreter)      = false
CC.may_discard_trees(interp::JETInterpreter) = false

# JETInterpreter specific
# -----------------------

JETAnalysisParams(interp::JETInterpreter) = interp.analysis_params

JETLogger(interp::JETInterpreter) = interp.logger

# TODO do report filtering or something configured by `JETAnalysisParams(interp)`
function report!(interp::JETInterpreter, report::InferenceErrorReport)
    push!(interp.reports, report)
end

function stash_uncaught_exception!(interp::JETInterpreter, report::UncaughtExceptionReport)
    push!(interp.uncaught_exceptions, report)
end

# check if we're in a toplevel module
@inline istoplevel(interp::JETInterpreter, sv::InferenceState)    = istoplevel(interp, sv.linfo)
@inline istoplevel(interp::JETInterpreter, linfo::MethodInstance) = interp.toplevelmod === linfo.def

@inline istoplevelmod(interp::JETInterpreter, mod::Module) = mod in interp.toplevelmods

is_global_slot(interp::JETInterpreter, slot::Int)   = slot in keys(interp.global_slots)
is_global_slot(interp::JETInterpreter, slot::Slot)  = is_global_slot(interp, slot_id(slot))
is_global_slot(interp::JETInterpreter, sym::Symbol) = sym in values(interp.global_slots)

@inline with_toplevel_logger(@nospecialize(f), interp::JETInterpreter, @nospecialize(filter = ≥(DEFAULT_LOGGER_LEVEL)); kwargs...) =
    with_logger(f, JETLogger(interp).toplevel_logger, filter, "toplevel"; kwargs...)

@inline with_inference_logger(@nospecialize(f), interp::JETInterpreter, @nospecialize(filter = ≥(DEFAULT_LOGGER_LEVEL)); kwargs...) =
    with_logger(f, JETLogger(interp).inference_logger, filter, "inference"; kwargs...)

@inline function with_logger(
    @nospecialize(f), io::Union{Nothing,IO}, @nospecialize(filter), logger_name;
    @nospecialize(pre = identity))
    isnothing(io) && return
    level = get_logger_level(io)
    filter(level) || return
    pre(io)
    print(io, "[$logger_name-$(LOGGER_LEVELS[level])] ")
    f(io)
end


text:
using BinaryProvider
import Libdl

version = v"0.2.9"

if arch(platform_key_abi()) != :x86_64
    @warn "FastTransforms has only been tested on x86_64 architectures."
end

const extension = Sys.isapple() ? "dylib" : Sys.islinux() ? "so" : Sys.iswindows() ? "dll" : ""

print_error() = error(
    "FastTransforms could not be properly installed.\n Please check that you have all dependencies installed. " *
    "Sample installation of dependencies:\n" *
    print_platform_error(platform_key_abi())
)

print_platform_error(p::Platform) = "On $(BinaryProvider.platform_name(p)), please consider opening a pull request to add support.\n"
print_platform_error(p::MacOS) = "On MacOS\n\tbrew install gcc@8 fftw mpfr\n"
print_platform_error(p::Linux) = "On Linux\n\tsudo apt-get install gcc-8 libblas-dev libopenblas-base libfftw3-dev libmpfr-dev\n"
print_platform_error(p::Windows) = "On Windows\n\tvcpkg install openblas:x64-windows fftw3[core,threads]:x64-windows mpir:x64-windows mpfr:x64-windows\n"

# Rationale is as follows: The build is pretty fast, so on Linux it is typically easiest
# to just use the gcc of the system to build the library and include it. On MacOS, however,
# we need to actually install a gcc first, because Apple's OS comes only shipped with clang,
# so here we download the binary.
ft_build_from_source = get(ENV, "FT_BUILD_FROM_SOURCE", Sys.isapple() ? "false" : "true")
if ft_build_from_source == "true"
    make = Sys.iswindows() ? "mingw32-make" : "make"
    compiler = Sys.isapple() ? "CC=gcc-8" : "CC=gcc"
    flags = Sys.isapple() ? "FT_USE_APPLEBLAS=1" : Sys.iswindows() ? "FT_FFTW_WITH_COMBINED_THREADS=1" : ""
    script = """
        set -e
        set -x
        if [ -d "FastTransforms" ]; then
            cd FastTransforms
            git fetch
            git checkout v$version
            cd ..
        else
            git clone -b v$version https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
        fi
        cd FastTransforms
        $make lib $compiler $flags
        cd ..
        mv -f FastTransforms/libfasttransforms.$extension libfasttransforms.$extension
    """
    try
        run(`bash -c $(script)`)
    catch
        print_error()
    end
    println("FastTransforms built from source.")
else
    const GCC = BinaryProvider.detect_compiler_abi().gcc_version
    namemap = Dict(:gcc4 => "gcc-4.9", :gcc5 => "gcc-5", :gcc6 => "gcc-6",
                   :gcc7 => "gcc-7", :gcc8 => "gcc-8", :gcc9 => "gcc-9")
    if !(GCC in keys(namemap))
        error("Please ensure you have a version of gcc from gcc-4.9 to gcc-9.")
    end
    try
        download("https://github.com/MikaelSlevinsky/FastTransforms/releases/download/" *
                 "v$version/libfasttransforms.v$version.$(namemap[GCC]).$extension",
                 joinpath(dirname(@__DIR__), "deps", "libfasttransforms.$extension"))
    catch
        print_error()
    end
    println("FastTransforms installed by downloading binaries.")
end