text:
function total_derivative(params_baseline, ϵ = 0.000001, settings = settings_defaults())
    @unpack ρ, d, γ = params_baseline
    sol_baseline = stationary_algebraic(params_baseline, settings)
    if γ == 1
        U_1 = 1/(ρ*sol_baseline.c)
        U_2 = (1/ρ^2)
    else
        U_1 = sol_baseline.c^(-γ) / (ρ + (γ-1)* sol_baseline.g)
        U_2 = sol_baseline.c^(1-γ) * (ρ + (γ-1)* sol_baseline.g)^(-2)
    end

    d_counterfactual = d + ϵ
    params_counterfactual = merge(params_baseline, (d = d_counterfactual,))
    sol_counterfactual = stationary_algebraic(params_counterfactual, settings)

	ACR = (-1/params_baseline.θ)*log(sol_counterfactual.λ_ii/sol_baseline.λ_ii)

    sol_fc_d = steady_state_from_g(sol_baseline.g, sol_baseline.z_hat, sol_baseline.Ω, params_counterfactual, settings)
    sol_fc_g = steady_state_from_g(sol_baseline.g + ϵ, sol_baseline.z_hat, sol_baseline.Ω, params_baseline, settings)
    sol_fc_zhat = steady_state_from_g(sol_baseline.g, sol_baseline.z_hat + ϵ, sol_baseline.Ω, params_baseline, settings)
    sol_fc_Omega = steady_state_from_g(sol_baseline.g, sol_baseline.z_hat, sol_baseline.Ω + ϵ, params_baseline, settings)

    partial_fc_d = (sol_fc_d.c - sol_baseline.c)/(d_counterfactual - params_baseline.d)
    partial_fc_g = (sol_fc_g.c - sol_baseline.c)/(sol_fc_g.g - sol_baseline.g)
    partial_fc_zhat = (sol_fc_zhat.c - sol_baseline.c)/(sol_fc_zhat.z_hat - sol_baseline.z_hat)
    partial_fc_Omega = (sol_fc_Omega.c - sol_baseline.c)/(sol_fc_Omega.Ω - sol_baseline.Ω)

    partial_Omega_d = (sol_counterfactual.Ω - sol_baseline.Ω)/(d_counterfactual - params_baseline.d)
    partial_zhat_d = (sol_counterfactual.z_hat - sol_baseline.z_hat)/(d_counterfactual - params_baseline.d)
    partial_g_d = (sol_counterfactual.g - sol_baseline.g)/(d_counterfactual - params_baseline.d)

    d_U_d_total = (sol_counterfactual.U_bar - sol_baseline.U_bar)/(d_counterfactual - params_baseline.d)
    total_decomp = U_1*(partial_fc_Omega*partial_Omega_d + partial_fc_zhat*partial_zhat_d + partial_fc_g*partial_g_d + partial_fc_d) + U_2*partial_g_d
    check = total_decomp/d_U_d_total
	total_decomp_CE = ρ*total_decomp

    decomp_fc_d = partial_fc_d
    decomp_fc_Omega_Omega_d = partial_fc_Omega*partial_Omega_d
    decomp_fc_zhat_zhat_d = partial_fc_zhat*partial_zhat_d
    decomp_fc_g_g_d = partial_fc_g*partial_g_d
    decomp_g_d = partial_g_d

    U1_partial_fc_d_frac = U_1*partial_fc_d/total_decomp
    U1_decomp_fc_Omega_Omega_d_frac = U_1*decomp_fc_Omega_Omega_d/total_decomp
	U1_decomp_fc_zhat_zhat_d_frac = U_1*decomp_fc_zhat_zhat_d/total_decomp
	U1_decomp_fc_g_g_d_frac = U_1*decomp_fc_g_g_d/total_decomp
	U2_decomp_g_d_frac = U_2*decomp_g_d/total_decomp

    planner_0_g_CE = ρ*(U_1*partial_fc_g +U_2)
    planner_0_Omega_CE = ρ*U_1*partial_fc_Omega
    planner_0_zhat_CE = ρ*U_1*partial_fc_zhat
	MU_g_CE = (1/ρ)
	semi_elasticity_g_d = params_baseline.d*partial_g_d
	U1_partial_fc_d_ACRfrac = U1_partial_fc_d_ACRfrac = ρ*U_1*partial_fc_d*ϵ/ACR
	d_T = 1 + 0.9*(params_baseline.d - 1)
	d_d = d_T - params_baseline.d
	CE_extrap10p = total_decomp_CE * d_d
	g_extrap10p = partial_g_d * d_d

    return (check = check, d_U_d_total = d_U_d_total, U_1 = U_1, U_2 = U_2,
										ACR = ACR,
										∂_fc_d = partial_fc_d,
										∂_fc_g = partial_fc_g,
										∂_fc_zhat = partial_fc_zhat,
										∂_fc_Omega = partial_fc_Omega,
										∂_g_d = partial_g_d,
										∂_zhat_d = partial_zhat_d,
										∂_Omega_d = partial_Omega_d,
										total_decomp = total_decomp,
										decomp_fc_d = decomp_fc_d,
										decomp_fc_Omega_Omega_d = decomp_fc_Omega_Omega_d ,
										decomp_fc_zhat_zhat_d = decomp_fc_zhat_zhat_d,
										decomp_fc_g_g_d = decomp_fc_g_g_d,
										decomp_g_d = decomp_g_d,
										U1_partial_fc_d_frac = U1_partial_fc_d_frac,
										U1_decomp_fc_Omega_Omega_d_frac = U1_decomp_fc_Omega_Omega_d_frac,
										U1_decomp_fc_zhat_zhat_d_frac = U1_decomp_fc_zhat_zhat_d_frac,
										U1_decomp_fc_g_g_d_frac = U1_decomp_fc_g_g_d_frac,
										U2_decomp_g_d_frac = U2_decomp_g_d_frac,
										planner_0_g_CE = planner_0_g_CE,
										planner_0_Omega_CE = planner_0_Omega_CE,
										planner_0_zhat_CE = planner_0_zhat_CE,
										MU_g_CE = MU_g_CE,
										semi_elasticity_g_d = semi_elasticity_g_d,
										U1_partial_fc_d_ACRfrac = U1_partial_fc_d_ACRfrac,
										d_T = d_T,
										d_d = d_d,
										CE_extrap10p = CE_extrap10p,
										g_extrap10p = g_extrap10p)
end

function steady_state_from_g(g, z_hat, Ω, parameters, settings)
    @assert parameters.υ > 0 && parameters.κ > 0
    @assert z_hat > 1 && Ω > 0 && g > 0 # input validation still required
    staticvalues = staticvals([g, z_hat, Ω], parameters)
    return merge(staticvalues, merge((g = g, z_hat = z_hat, Ω = Ω,), welfare([g, z_hat, Ω], parameters)))
end

function stationary_algebraic(parameters, settings)
    x0 = settings.stationary_x0(parameters, settings)
    @assert parameters.υ > 0 && parameters.κ > 0
    sol = nlsolve(x -> stationary_algebraic_aux(x, parameters), x0, inplace = false) # uses distorted parameters
    converged(sol) || error("Solver didn't converge.")
    g, z_hat, Ω = sol.zero
    @assert z_hat > 1 && Ω > 0 && g > 0
    staticvalues = staticvals([g, z_hat, Ω], parameters) # uses undistorted parameters (except x)

    # calculate value function using differential objects based on the grid
    @unpack z_ex = settings
    z = z_ex[2:end-1]
    @unpack θ, σ = parameters
    bc = (Mixed(ξ = σ-1), Mixed(ξ = σ-1)) # boundary conditions for differential operators
    L_1_minus = L₁₋bc(z_ex, bc) # use backward difference as the drift is negative
    L_2 = L₂bc(z_ex, bc)
    ω = ω_weights(z_ex, θ, σ-1) # get quadrature weights for the distribution on the rescaled grid.
    Ξ₁ = 1/(1 - (σ-1)*(z[1] - z_ex[1])) # (24), with ξ = (σ-1)
    p = (L_1_minus = L_1_minus, L_2 = L_2, ω = ω, Ξ₁ = Ξ₁, z = z)
    v_tilde = stationary_numerical_given_vals([g, z_hat, Ω], p, parameters, settings).v_tilde

    return merge(staticvalues, merge((g = g, z_hat = z_hat, Ω = Ω, v_tilde = v_tilde), welfare([g, z_hat, Ω], parameters))) # calculate quantities and return
end

function welfare(vals, parameters)
    g, z_hat, Ω = vals
    @unpack ρ, σ, N, θ, γ, d, κ, ζ, η, Theta, χ, υ, μ, δ = parameters
    @unpack F, r, ν, a, b, S, L_tilde, z_bar, w, π_min = staticvals(vals, parameters)
    c_T = c(L_tilde, Ω, z_bar)
    return (y = c_T, # y = c by (C.78)
            c = c_T,
            λ_ii = λ_ii(z_hat, parameters),
            U_bar = γ == 1 ? (ρ*log(c_T) + g) / ρ^2 : 1/(1-γ) * (c_T^(1-γ))/(ρ + (γ-1)*g)) # (C.16)
end

# Kernel for the algebraic stationary.
function stationary_algebraic_aux(vals, parameters)
    @unpack ρ, σ, N, θ, γ, d, κ, η, Theta, υ, μ, δ, χ, ζ = parameters
    @unpack F, r, ν, a, b, S, L_tilde, z_bar, w, x, π_min = staticvals(vals, parameters)
    g, z_hat, Ω = vals

    big_denom = ν*(θ + ν)*(θ - σ + 1) # (C.19)
    denom_1 = a*(g - r) # (C.19)
    num_1 = ν*(N-1)*(θ - σ + 1)*(d^(1 - σ)*(θ + ν)*z_hat^(-θ + σ - 1)-b*θ*z_hat^(-θ-ν)) # (C.19)
    num_2 = θ*(ν*(N-1)*d^(1-σ)*(θ+ν)*z_hat^(-θ + σ -1) + (ν + σ - 1)*(θ + ν - σ + 1)) # (C.19)
    return [x/π_min - a*(χ/(1-χ))*(σ + ν - 1)/ν, # (C.18)
            1 + (σ-1)/ν - (num_1/denom_1 + num_2)/big_denom + (χ/(1-χ))*(σ + ν - 1)/(ν), # (C.19)
            π_min - (1- L_tilde)/((σ -1)*z_bar^(σ-1))] # (C.20)
end

function staticvals(vals, parameters)
    g, z_hat, Ω = vals
    @assert parameters.η ≈ 0  # we think that the L_tilde formulas are wrong for η != 0
    @unpack ρ, σ, N, θ, γ, d, κ, ζ, η, Theta, χ, υ, μ, δ = parameters

    F(z) = 1 - z^(-θ) # (C.1)
    r = ρ + γ*g + δ # (C.6)
    ν = (μ-g)/υ^2 + sqrt(((g-μ)/υ^2)^2 + (r-g)/(υ^2/2)) # (C.3)
    a = (r - g - (σ - 1)*(μ - g + (σ - 1)*υ^2/2))^(-1) # (C.4)
    b = (1 - a*(r-g))*d^(1-σ)*z_hat^(ν + σ - 1) # (C.5)
    S_T = S(g, parameters) # calls like this rely on utils/static.jl
    L_tilde_T = L_tilde(g, z_hat, Ω, δ, S_T, parameters)
    L_tilde_x_T = L_tilde_x(z_hat, Ω, parameters)
    L_tilde_E_T = L_tilde_E(Ω, δ, parameters)
    L_tilde_a_T = L_tilde_a(Ω, S_T, parameters)
    z_bar_T = z_bar(z_hat, Ω, parameters)
    w_T = w(z_bar_T, parameters)
    x_T = x(w_T, parameters)
    π_min_T = (d^(σ-1) * κ)/(z_hat^(σ-1)) # (C.12, inverted to express π_min as a function of parameters and z_hat) This is where the export threshold is used.  pi_min(L_tilde, z_bar)  is already used in the system of equations
    π_rat_T = π_rat(z_hat, parameters)

    return (F = F, r = r, ν = ν, a = a, b = b, S = S_T, L_tilde = L_tilde_T, L_tilde_x = L_tilde_x_T, L_tilde_E = L_tilde_E_T, L_tilde_a = L_tilde_a_T,
            z_bar = z_bar_T, w = w_T, x = x_T, π_min = π_min_T, π_rat = π_rat_T)
end

# Kernel of the stationary numerical.
function stationary_numerical_given_vals(vals, p, parameters, settings)
    g, z_hat, Ω = vals
    @unpack F, r, ν, a, b, S, L_tilde, z_bar, w, π_min = staticvals([g, z_hat, Ω], parameters)
    @unpack ρ, σ, N, θ, γ, d, κ, η, Theta, ζ, χ, υ, μ, δ = parameters
    @unpack Ξ₁, L_1_minus, L_2, ω, z = p
    r_tilde = r - g - 0 # (C.59, and g_w = 0 at steady state)
    ρ_tilde = r_tilde - (σ - 1)*(μ - g + (σ-1)*(υ^2/2)) # (C.41)
    A_T = ρ_tilde * I - (μ - g + (σ-1)*υ^2)*L_1_minus - υ^2/2 * L_2 # (52)
    i(z) = z >= log(z_hat) ? 1 : 0 # indicator function for next equation.
    π(z) = π_min * (1 + (N-1)*d^(1-σ)*i(z)) - (N-1)*κ*exp(-(σ-1)*z)*i(z) # (39)
    v_tilde = A_T \ π.(z) # discretized system of ODE for v, where v'(T) = 0 (53)

    # System of equations, given the numerically solved ODE
    return (residuals = [Ξ₁*v_tilde[1] - dot(v_tilde, ω) + ζ, # (54)
            Ξ₁*v_tilde[1] - ζ*(1-χ)/χ, # (56)
            π_min - (1 - L_tilde)/((σ-1)*z_bar^(σ-1))], # (38)
            v_tilde = v_tilde)
end

# Solve the ODE by root-finding. Should use the analytic most of the time, but this is a useful check.
function stationary_numerical(parameters, settings)
    @unpack z_ex = settings
    x0 = settings.stationary_x0(parameters, settings)
    @unpack ρ, σ, N, θ, γ, d, κ, η, Theta, υ, μ, δ = parameters
    @assert υ > 0 && κ > 0 # Parameter validation
    z = z_ex[2:end-1] # form a uniform extended grid

    # Differential objects we need for the residuals
    bc = (Mixed(ξ = σ-1), Mixed(ξ = σ-1)) # boundary conditions for differential operators
    L_1_minus = L₁₋bc(z_ex, bc) # use backward difference as the drift is negative
    L_2 = L₂bc(z_ex, bc)
    ω = ω_weights(z_ex, θ, σ-1) # Get quadrature weights for the distribution on the rescaled grid.
    Ξ₁ = 1/(1 - (σ-1)*(z[1] - z_ex[1])) # (24), with ξ = (σ-1)

    p = (Ξ₁ = Ξ₁, L_1_minus = L_1_minus, L_2 = L_2, ω = ω, z = z)
    sol = nlsolve(vals -> stationary_numerical_given_vals(vals, p, parameters, settings).residuals, x0, inplace = false)
    converged(sol) || error("Solver didn't converge.")
    g_T, z_hat_T, Ω_T = sol.zero

    staticvalues = staticvals([g_T, z_hat_T, Ω_T], parameters)
    v_tilde = stationary_numerical_given_vals([g_T, z_hat_T, Ω_T], p, parameters, settings).v_tilde
    return merge(staticvalues, merge((g = g_T, z_hat = z_hat_T, Ω = Ω_T, v_tilde = v_tilde), welfare([g_T; z_hat_T; Ω_T], parameters)))
end

# Utilities for comparison, display, etc.
function compare_steady_states(parameters, settings; verbose = false, algebraic = false)

    if algebraic
        stationary = stationary_algebraic
    else
        stationary = stationary_numerical
    end

    parameters_0 = merge(parameters, (d = parameters.d_0,))
    parameters_T = merge(parameters, (d = parameters.d_T,))
    stationary_0 = stationary(parameters_0, settings)
    stationary_T = stationary(parameters_T, settings)
    change_welfare = 100*consumption_equivalent(stationary_T.U_bar, stationary_0.U_bar, parameters)
    change_trade = (1-stationary_T.λ_ii) - (1-stationary_0.λ_ii)

    if verbose
        println("SS to SS welfare gain: $(change_welfare)")
        println("Change in Trade: $(change_trade)")
        println("Growth Rates across SS: $(100*stationary_T.g) (time T) vs $(100*stationary_0.g) (time 0)")
    end

    return (stationary_0 = stationary_0,
            stationary_T = stationary_T,
            change_welfare = change_welfare,
            change_trade = change_trade)
end


text:
"""
    natural_eigdist(𝚽, 𝛌, Q; α = 1.0, T = :Inf,
                    input_format = :zero_measures, distance = :DAG,
                    edge_weight = 1, edge_length = 1)

compute natural distances between graph Laplacian eigenvectors.

# Input Arguments
- `𝚽::Matrix{Float64}`: matrix of (weighted) graph Laplacian eigenvectors.
- `𝛌::Vector{Float64}`: vector of eigenvalues.
- `Q::Matrix{Float64}`: unweighted incidence matrix of the graph.
- `α::Float64`: ROT parameter. (default: `1.0`)
- `T::Any`: TSD parameter, i.e., the stopping time T in K_functional (default: `:Inf`)
- `input_format::Symbol`: options: `:zero_measures`, `:pmf1` and `:pmf2` (default: `:zero_measures`)
- `distance::Symbol`: options: `:ROT`, `:HAD`, `:DAG` and `:TSD` (default: `:DAG`)
- `edg_length::Any`: vector of edge lengths (default: 1 represents unweighted graphs)
- `edge_weight::Any`: the weights vector, which stores the affinity weight of
    each edge (default: `1` represents unweighted graphs).

# Output Argument
- `dis::Matrix{Float64}`: the distance matrix, dis[i,j] = d(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

"""
function natural_eigdist(𝚽, 𝛌, Q; α = 1.0, T = :Inf,
                         input_format = :zero_measures, distance = :DAG,
                         edge_weight = 1, edge_length = 1)
    N = size(Q, 1)
    P = deepcopy(𝚽)
    if input_format == :zero_measures
        P[:, 1] .= 0
    elseif input_format == :pmf1
        P = P.^2
    elseif input_format == :pmf2
        P = exp.(P) ./ sum(exp.(P), dims = 1)
    else
        @error("input_format does not support $(input_format)!")
        return
    end

    if distance == :ROT
        D = eigROT_Distance(P, Q; edge_length = edge_length, α = α)
    elseif distance == :HAD
        D = eigHAD_Distance(𝚽, 𝛌)
    elseif distance == :DAG
        D = eigDAG_Distance(𝚽, Q, N; edge_weight = edge_weight)
    elseif distance == :TSD
        tL = Q * Q'
        t𝛌, t𝚽 = eigen(Matrix(tL))
        D = eigTSD_Distance(P, t𝚽, t𝛌, Q; length = edge_length, T = T)
    else
        error("distance does not support $(distance)!")
        return
    end

    return D
end

ext: jl


text:
# Use this version of Elemental
Elsha = "79987d38b04838acf6b6195be1967177521ee908"

using Compat

if is_windows()
    error("Elemental only works on Unix Platforms")
end

depdir = dirname(@__FILE__)

if !isdir(joinpath(depdir, "src"))
    mkdir(joinpath(depdir, "src"))
end
srcdir = joinpath(depdir, "src", "Elemental")

if !isdir(joinpath(depdir, "usr"))
    mkdir(joinpath(depdir, "usr"))
end
prefix = joinpath(depdir, "usr")

if VERSION < v"0.5.0-dev+5398"
    if !isdir(srcdir)
        Base.Git.run(`clone -- https://github.com/elemental/Elemental.git $srcdir`)
    end
    cd(srcdir) do
        Base.Git.run(`checkout $Elsha`)
    end
else
    if !isdir(srcdir)
        LibGit2.clone("https://github.com/elemental/Elemental.git", "$srcdir")
    end
    cd(srcdir) do
        LibGit2.checkout!(LibGit2.GitRepo("."), "$Elsha")
    end
end

if VERSION < v"0.5.0-dev+4343"
    Base.check_blas()
    blas = Base.blas_vendor()
else
    BLAS.check()
    blas = BLAS.vendor()
end
mathlib = Libdl.dlpath(BLAS.libblas)
blas64 = LinAlg.USE_BLAS64 ? "ON" : "OFF"
blas_suffix = blas === :openblas64 ? "_64_" : "_"
build_procs = (haskey(ENV, "CI") && ENV["CI"] == "true") ? 2 : Sys.CPU_CORES

builddir = joinpath(depdir, "builds")
if isdir(builddir)
    rm(builddir, recursive=true)
end
mkdir(builddir)

cd(builddir) do
    run(`cmake -D CMAKE_INSTALL_PREFIX=$prefix
               -D INSTALL_PYTHON_PACKAGE=OFF
               -D PYTHON_EXECUTABLE=""
               -D PYTHON_SITE_PACKAGES=""
               -D EL_USE_64BIT_INTS=$blas64
               -D EL_USE_64BIT_BLAS_INTS=$blas64
               -D MATH_LIBS=$mathlib
               -D EL_BLAS_SUFFIX=$blas_suffix
               -D EL_LAPACK_SUFFIX=$blas_suffix
               -D CMAKE_INSTALL_RPATH=$prefix/lib
               $srcdir`)
    run(`make -j $build_procs`)
    run(`make install`)
end


text:
# module for testing each slam modules

module TestSlam
  using Test

  # target modules
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/estimated_object.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/map_particle.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/fast_slam_1/anime_fast_slam_1.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/fast_slam_2/anime_fast_slam_2.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/graph_based_slam/psi_sensor.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/graph_based_slam/logger_agent.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/graph_based_slam/traj_obsrv_logger.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/slam/graph_based_slam/anime_graph_based_slam.jl"))

  function main()
    @testset "Slam" begin
      @testset "EstimatedObject" begin
        ep = EstimatedObject(1.0, 2.0, shape=:circle, color=:red,
                             size=5, id=20, cov=[1 1;2 2])
        @test ep.pose == [1.0, 2.0]
        @test ep.id == 20
        @test ep.shape == :circle
        @test ep.color == :red
        @test ep.size == 5
        @test ep.cov == [1 1;2 2]
      end
      @testset "MapParticle" begin
        mp = MapParticle([1.0 2.0 3.0], 0.4, 30)
        @test mp.pose == [1.0 2.0 3.0]
        @test mp.weight == 0.4
        @test mp.object_num == 30
        @test length(mp.map.objects) == 30
      end
      @testset "FastSlam1.0" begin
        @test_nowarn AnimeFastSlam1.main(is_test=true)
      end
      @testset "FastSlam2.0" begin
        @test_nowarn AnimeFastSlam2.main(is_test=true)
      end
      @testset "PsiSensor" begin
        m = Map()
        add_object(m, Object(2.0, -2.0))
        s = PsiSensor(m)
        @test s.last_data == []
        @test s.dist_rng == (0.5, 6.0)
        @test s.dir_rng == (-pi/3, pi/3)
        @test visible(s, [0.5, -pi/3]) == true
        @test visible(s, [6.0, pi/3]) == true
        @test visible(s, [0.2, -pi/3]) == false
        @test visible(s, [0.5, pi]) == false
      end
      @testset "LoggerAgent" begin
        a = LoggerAgent(0.1, 1.0, is_test=true)
        obsrv = []
        @test a.pose == [0.0, 0.0, 0.0]
        speed, yaw_rate = draw_decision!(a, obsrv)
        @test speed == 0.1
        @test yaw_rate == 1.0
        @test a.time_interval == 0.1
        @test a.estimator === nothing
        @test a.prev_spd == 0.0
        @test a.prev_yr == 0.0
        @test a.step == 0
      end
      @testset "TrajObsrvLogger" begin
        @test_nowarn TrajObsrvLogger.main(is_test=true)
      end
      @testset "GraphBasedSlam" begin
        @test_nowarn AnimeGraphBasedSlam.main(true)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
Base.iterate(tree::Tree) = iterate((tree, x->true, typemax(Int)))

function Base.iterate(tuple::Tuple{Tree, Function, Int, Int}, state = (tuple[1], 0, Vector{Int}()))
    cell, count, indices = state

    if cell == nothing return nothing end
    tree, filter, min_level, max_level = tuple

    cell, next_cell = find_next_filtered_cell!(cell, indices, filter, min_level, max_level)
    if cell == nothing
        return nothing
    else
        return (cell, (next_cell, count + 1, indices))
    end

end

function find_next_filtered_cell!(cell::Tree{N}, indices::Vector{Int}, filter::Function, min_level::Int, max_level::Int) where N
    next_cell = find_next_cell!(cell, indices, max_level)
    while !filter(cell) || (min_level > 0 && level(cell) < min_level)
        cell = next_cell
        if cell == nothing break end
        next_cell = find_next_cell!(cell, indices, max_level)
    end

    return cell, next_cell
end

# indices[l] yields the child index at level l in the following sense
# cell = cell.parent.children[indices[end]]
#         = cell.parent.parent.children[indices[end-1]].children[end]
#         = ...
function find_next_cell!(cell::Tree{N}, indices::Vector{Int}, max_level::Int) where N
    if !active(cell) && level(cell) < max_level
        # Depth first: go to higher level
        next_cell = cell.children[1]
        push!(indices, 1)
    elseif isempty(indices)
        next_cell = nothing
    else
        # Find sibling or sibling of parent (or grandparent etc.)
        next_cell = cell
        while indices[end] == 1<<N
            next_cell = next_cell.parent
            pop!(indices)
            if isempty(indices) return nothing end
        end
        indices[end] += 1
        next_cell = next_cell.parent.children[indices[end]]
    end

    return next_cell
end

function Base.show(io::IO, tree::Tree)
    compact = get(io, :compact, false)

    print(io, "$(typeof(tree)) ")
    if initialized(tree)
        if level(tree) == 0
            print(io, "root ")
        else
            if active(tree) print(io, "leaf ") end
            print(io, "on level $(level(tree)) ")
        end
        if !active(tree) && !compact
            print(io, "with $(levels(tree)) levels and $(length(cells(tree, filter=active))) active_cells out of $(length(cells(tree))) cells")
        end
    end
end

@inline cells(tree::Tree; filter::Function = cell -> true, min_level = 0, max_level = typemax(Int)) = (tree, filter, min_level, max_level)
@inline cells(tree::Tree, level::Int; filter::Function = cell -> true) = (tree, filter, level, level)

function Base.length(tuple::Tuple{Tree, Function, Int, Int})
    count = 0
    for cell ∈ tuple
        count += 1
    end
    return count
end

# Iterate over the faces, starting at the left-hand side (first dimension) face of some tree
function Base.iterate(tuple::Tuple{Face, Function, Int, Int}, state = (tuple[1], 0, tuple[1].cells[2], Vector{Int}(), [1, 1]))
    face, count, cell, cell_indices, face_indices = state

    if face == nothing return nothing end
    tree, filter, min_level, max_level = tuple

    next_face, cell = find_next_face!(face_indices, cell_indices, cell, max_level)
    face_good = filter(face) && (min_level == 0 || level(face) >= min_level)
    while cell != nothing && !face_good
        face = next_face
        next_face, cell = find_next_face!(face_indices, cell_indices, cell, max_level)
        face_good = filter(face) && (min_level == 0 || level(face) >= min_level)
    end
    if cell == nothing next_face = nothing end

    if !face_good
        return nothing
    else
        return (face, (next_face, count + 1, cell, cell_indices, face_indices))
    end
end

# Loop over the faces (NB now the filter acts on faces)
# Per cell we loop over the left-hand side faces as well as boundary and refinement faces
function find_next_face!(face_indices::Vector{Int}, cell_indices::Vector{Int}, cell::Tree{N}, max_level::Int) where N
    if face_indices[2] == 1
        # Go to right hand side
        face_indices[2] += 1
        if regular(cell.faces[face_indices[1],face_indices[2]])
            # but only if face is not regular, otherwise go to next dimension
            face_indices[2] = 1
            face_indices[1] += 1
        end
    else
        # Go to next dimension
        face_indices[2] = 1
        face_indices[1] += 1
    end
    if face_indices[1] > N
        cell = find_next_cell!(cell, cell_indices, max_level)
        if cell == nothing return nothing, nothing end
        face_indices .= 1
    end
    return cell.faces[face_indices[1],face_indices[2]], cell
end

# Iterators over all of the faces
@inline all_faces(tree::Tree; filter::Function = face -> true, min_level::Int = 0, max_level::Int = typemax(Int)) = (tree.faces[1], filter, min_level, max_level)
@inline all_faces(tree::Tree, level::Int; filter::Function = face -> true) = (tree.faces[1], filter, level, level)
function Base.length(tuple::Tuple{Face, Function, Int, Int})
    count = 0
    for cell ∈ tuple
        count += 1
    end
    return count
end


text:
@testset "Events/CDATA Marked Section" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    CDataStart = E.CDATAMarkedSectionStart
    CDataEnd   = E.CDATAMarkedSectionEnd

    DC = E.DataContent
    ME = E.MarkupError

    @testset "Events/CDATA Marked Section (Positive)" begin
        @test evaluate("<![CDATA[]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                            DC("", L.Location("a buffer", -1)),
                                            CDataEnd(L.Location("a buffer", -1)) ]

        @test evaluate("<![CDATA[a]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                             DC("a", L.Location("a buffer", -1)),
                                             CDataEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![CDATA[abc]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                               DC("abc", L.Location("a buffer", -1)),
                                               CDataEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![CDATA[é]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                             DC("é", L.Location("a buffer", -1)),
                                             CDataEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![CDATA[aé]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                              DC("aé", L.Location("a buffer", -1)),
                                              CDataEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![CDATA[éb]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                              DC("éb", L.Location("a buffer", -1)),
                                              CDataEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![CDATA[aéb]]>") == [ CDataStart(L.Location("a buffer", -1)),
                                               DC("aéb", L.Location("a buffer", -1)),
                                               CDataEnd(L.Location("a buffer", -1)) ]

        # This is from the XML specification, § 2.7.
        #
        @test (evaluate("<![CDATA[<greeting>Hello, world!</greeting>]]>")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("<greeting>Hello, world!</greeting>", L.Location("a buffer", -1)),
                    CDataEnd(L.Location("a buffer", -1)) ])
    end

    @testset "Events/CDATA Marked Section (Negative ... no section type (CDATA))" begin
        @test (evaluate("<![")
               == [ ME("ERROR: Expecting the start of a CDATA or conditional marked section.", L.Location("a buffer", -1)) ])

        # White space is not allowed between the DSO and the CDATA token.
        #
        @test (evaluate("<![ CDATA")
               == [ ME("ERROR: White space is not allowed before the 'CDATA' keyword.", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '[' to open a CDATA marked section.", L.Location("a buffer", -1)) ])

        @test (evaluate("<![ ")
               == [ ME("ERROR: Expecting the start of a CDATA or conditional marked section.", L.Location("a buffer", -1)),
                    DC(" ", true, L.Location("a buffer", -1)) ])

        # Check that a random token is caught.
        #
        events = evaluate("<![<")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting the start of a CDATA or conditional marked section.",
                              L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/CDATA Marked Section (Negative ... wrong section type (CDATA))" begin
        @test (evaluate("<![SDATA")
               == [ ME("ERROR: Expecting the start of a CDATA or conditional marked section.", L.Location("a buffer", -1)),
                    DC("SDATA", false, L.Location("a buffer", -1)) ])

        @test (evaluate("<![cdata")
               == [ ME("ERROR: The keyword 'cdata' must be uppercased.", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '[' to open a CDATA marked section.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/CDATA Marked Section (Negative ... no second DSO)" begin
        @test (evaluate("<![CDATA")
               == [ ME("ERROR: Expecting '[' to open a CDATA marked section.", L.Location("a buffer", -1)), ])

        # White space is not allowed between the CDATA token and the second DSO.
        #
        @test (evaluate("<![CDATA [")
               == [ ME("ERROR: White space is not allowed after the 'CDATA' keyword.", L.Location("a buffer", -1)),
                    CDataStart(L.Location("a buffer", -1)),
                    DC("", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting ']]>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)) ])
    end

    @testset "Events/CDATA Marked Section (Negative ... no MSC)" begin
        @test (evaluate("<![CDATA[")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting ']]>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)) ])

        @test (evaluate("<![CDATA[Hello, World!")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("Hello, World!", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting ']]>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)) ])
    end

    @testset "Events/CDATA Marked Section (Negative ... no TAGC)" begin
        @test (evaluate("<![CDATA[]]")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)) ])

        @test (evaluate("<![CDATA[Hello, World!]]")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("Hello, World!", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)) ])

        # White space is not allowed between the MSC adn the TAGC.
        #
        @test (evaluate("<![CDATA[]] >")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)),
                    DC(" ", true, L.Location("a buffer", -1)),
                    DC(">", false, L.Location("a buffer", -1)) ])

        @test (evaluate("<![CDATA[Hello, World!]] >")
               == [ CDataStart(L.Location("a buffer", -1)),
                    DC("Hello, World!", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end a CDATA marked section.", L.Location("a buffer", -1)),
                    CDataEnd(true, L.Location("a buffer", -1)),
                    DC(" ", true, L.Location("a buffer", -1)),
                    DC(">", false, L.Location("a buffer", -1)) ])
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
include("../src/BlueMesh7.jl")
using .BlueMesh7

using Random
using DataFrames
using Statistics: mean
using ProgressMeter
using BSON: @load, @save
using Flux
using CUDA
# ■

benchmark = DataFrame()
minutes = 7
ngraphs = 4

nactors = [32, 64, 96]
dims = [(25, 25), (30, 30), (35, 35)]

setups = [(d, n) for d in dims, n in nactors] |> vec
results = DataFrame()

Random.seed!(7)
allgraphs = [generate_positions(; n, dims) for _ = 1:ngraphs, (dims, n) in setups]

function bench(name, benchmark; init, unroll)
    println("Benching $(name)")
    bench_stats = DataFrame()

    @showprogress for (setup, graphs) in zip(setups, eachcol(allgraphs))
        setup_stats = DataFrame()
        spin = Threads.SpinLock()

        Threads.@threads for graph in graphs
            # GIL avoidance for python algorithms
            lock(spin)
            mesh = init(graph)
            unlock(spin)

            Random.seed!(7)
            stats = unroll(mesh, minutes)

            # since when it's not thread safe
            lock(spin)
            push!(setup_stats, stats)
            unlock(spin)
        end

        @assert nrow(setup_stats) == ngraphs
        averages = [setup, colwise(mean, setup_stats)...]
        push!(bench_stats, NamedTuple{(:setup, :pdr, :worstpdr, :delay, :centrality, :packetloss)}(averages))
    end

    bench_stats.algorithm = name

    for row in eachrow(bench_stats)
        push!(benchmark, row)
    end
end

bench("All relays (0)", benchmark,
      init = coordinates -> initialize_mesh(coordinates, ones(Int, length(coordinates))),
      unroll = (mesh, minutes) -> start(mesh; minutes))

bench("Half relays (0)", benchmark,
      init = coordinates -> initialize_mesh(coordinates, rand(0:1, length(coordinates))),
      unroll = (mesh, minutes) -> start(mesh; minutes))

# using PyCall
# pushfirst!(PyVector(pyimport("sys")."path"), "../test/ble-mesh-algorithms/standard-algorithms/")
# greedy = pyimport("greedy_algorithms")
# dominator = pyimport("dominator_algorithm")

# function init_rssi_to_adjacency(coordinates::Vector{NTuple{2, Int}}, Oracle)
#     protomesh = initialize_mesh(coordinates, ones(Int, length(coordinates)))
#     mask = protomesh.rssi_map .> protomesh.scanner_sensitivity + 1e-10

#     adjacency = zeros(Int, size(protomesh.rssi_map))
#     adjacency[mask] .= 1
#     adjacency = adjacency[2:end, 2:end]

#     roles = Oracle(adjacency)
#     initialize_mesh(coordinates, roles)
# end

# bench("Connect (1)", benchmark,
#       init = coordinates -> init_rssi_to_adjacency(coordinates, greedy.greedy_connect),
#       unroll = (mesh, minutes) -> start(mesh; minutes))

# bench("Dominator (1)", benchmark,
#       init = coordinates -> init_rssi_to_adjacency(coordinates, dominator.dominator),
#       unroll = (mesh, minutes) -> start(mesh; minutes))

@load "model.bson" model

function unroll_mc(env, minutes)
    moves = zeros(env.model.n)
    S = get_state(env)

    ps = model(S) |> softmax

    moves = rand.(mapslices(Categorical, ps', dims=2)) |> vec

    env(moves)

    for t = 1:minutes * 60_000
        env(moves, true)
    end

    getstats(env.model)
end

bench("Monte Carlo (3)", benchmark,
      init = (positions) -> BlueMesh7Env(positions),
      unroll = unroll_mc)

# ■

using Dates
datenow = Date(now())
@save "benchmark-$(datenow).bson" benchmark


text:
Infras = :(
module Infras
using MLStyle.MatchCore
using MLStyle.Extension
using MLStyle.Err
using MLStyle.Toolz: ($)
using MLStyle.Render: render, format

const __L__ = @__LINE__
export @format
macro format(args, template)
    args = Expr(:vect, :__L__ => LineNumberNode(__L__), :failed => QuoteNode(:($MatchCore.failed)), args.args...)
    esc(format(args, template))
end

export @typed_as
"""
this macro provide a convenient and type-stable way to generate
type checking.

However, the generated code requires a given variable `TARGET`
in metaprogramming level.

e.g

TARGET = :target_id
(@typed_as Int) # which generates a pattern to check if type is Int.

"""
macro typed_as(t)
    esc $ quote
        NAME = mangle(mod)
        __T__ = $t
        function (body)
            @format [body, tag, TARGET, NAME, __T__] quote

                @inline __L__ function NAME(TARGET :: __T__)
                    __T__
                    body
                end

                @inline __L__ function NAME(TARGET)
                    failed
                end

                NAME(tag)
            end
        end
    end
end

export patternOr
patternOr  = (p1, p2) -> body ->
    let p1 = p1(body), p2 = p2(body)
        tmp = mangle(Infras)
        @format [tmp, p1, p2] quote
            tmp = p1
            tmp === failed ? p2 : tmp
        end
    end

const APP_DESTRUCTORS = Vector{Tuple{Module, PDesc}}()
const GAPP_DESTRUCTORS = Vector{Tuple{Module, PDesc}}()

export def_pattern
function def_pattern(mod; predicate, rewrite, qualifiers=nothing)
    qualifiers = qualifiers === nothing ? Set([invasive]) : qualifiers
    desc = PDesc(predicate, rewrite, qualifiers)
    register_pattern(desc, mod)
end

# ======= App Patterns =============

"""
 To test if an datatype can
"""
function deconstructable
end

export register_app_pattern
function register_app_pattern(pdesc :: PDesc, def_mod::Module)
    push!(APP_DESTRUCTORS, (def_mod, pdesc))
end

export def_app_pattern
function def_app_pattern(mod; predicate, rewrite, qualifiers=nothing)
    qualifiers = qualifiers === nothing ? Set([invasive]) : qualifiers
    desc = PDesc(predicate, rewrite, qualifiers)
    register_app_pattern(desc, mod)
end

export mk_app_pattern
function mk_app_pattern(tag, hd, tl, use_mod)
    if isdefined(use_mod, hd)
        hd = getfield(use_mod, hd)
        for (def_mod, desc) in APP_DESTRUCTORS
            if qualifier_test(desc.qualifiers, use_mod, def_mod) && desc.predicate(hd, tl)
                return desc.rewrite(tag, hd, tl, use_mod)
            end
        end
    end
    info = string(hd) * "(" * join(map(string, tl), ", ") * ")"
    throw $ PatternUnsolvedException("invalid usage or unknown application case $info.")
end

# ===================================

# ===== Generalized App Patterns ====

export def_gapp_pattern
function def_gapp_pattern(mod; predicate, rewrite, qualifiers=nothing)
    qualifiers = qualifiers === nothing ? Set([invasive]) : qualifiers
    desc = PDesc(predicate, rewrite, qualifiers)
    register_gapp_pattern(desc, mod)
end

export register_gapp_pattern
function register_gapp_pattern(pdesc :: PDesc, def_mod::Module)
    push!(GAPP_DESTRUCTORS, (def_mod, pdesc))
end

export mk_gapp_pattern
function mk_gapp_pattern
end

# ===================================

def_pattern(Infras,
    predicate = x -> x isa Expr && x.head == :call,
    rewrite = (tag, case, mod) ->
    let hd = case.args[1], tl = case.args[2:end]
        @match hd begin
            ::Symbol   => mk_app_pattern(tag, hd, tl, mod)
            :($(::Symbol){$(_...)}) => mk_gapp_pattern(tag, [], hd, tl, mod)
            _ => throw("Deconstructor cannot be an expression like $(string(hd)).")
        end
    end)
end
)


text:
module SplittingPool

# scan into a sorted list of free buffers, splitting buffers along the way

using ..CuArrays
using ..CuArrays: @pool_timeit, @safe_lock

using DataStructures

using CUDAdrv

using Base: @lock
using Base.Threads: SpinLock

## tunables

# pool boundaries
const SMALL_CUTOFF   = 2^20    # 1 MiB
const LARGE_CUTOFF   = 2^25    # 32 MiB

# memory size rounding
const SMALL_ROUNDOFF = 2^9     # 512 bytes
const LARGE_ROUNDOFF = 2^17    # 128 KiB
const HUGE_ROUNDOFF  = 2^20    # 1 MiB

# maximum overhead (unused space) when selecting a buffer
# small and large buffers with large overhead will be split, but huge buffers never are
const SMALL_OVERHEAD = typemax(Int)
const LARGE_OVERHEAD = typemax(Int)
const HUGE_OVERHEAD  = 0

## block of memory

using Printf

@enum BlockState begin
    INVALID
    AVAILABLE
    ALLOCATED
    FREED
end

const block_id = Threads.Atomic{UInt}(1)

# TODO: it would be nice if this could be immutable, since that's what SortedSet requires
mutable struct Block
    ptr::Union{Nothing,CuPtr{Nothing}}  # base allocation
    sz::Int                             # size into it
    off::Int                            # offset into it

    state::BlockState
    prev::Union{Nothing,Block}
    next::Union{Nothing,Block}

    id::UInt

    Block(ptr, sz; off=0, state=INVALID, prev=nothing, next=nothing,
          id=Threads.atomic_add!(block_id, UInt(1))) =
        new(ptr, sz, off, state, prev, next, id)
end

Base.sizeof(block::Block) = block.sz
Base.pointer(block::Block) = block.ptr + block.off

iswhole(block::Block) = block.prev === nothing && block.next === nothing

## block utilities

function Base.show(io::IO, block::Block)
    fields = [@sprintf("#%d", block.id)]
    push!(fields, @sprintf("%s at %p", Base.format_bytes(sizeof(block)), pointer(block)))
    push!(fields, "$(block.state)")
    block.prev !== nothing && push!(fields, @sprintf("prev=Block(#%d)", block.prev.id))
    block.next !== nothing && push!(fields, @sprintf("next=Block(#%d)", block.next.id))

    print(io, "Block(", join(fields, ", "), ")")
end

# split a block at size `sz`, returning the newly created block
function split!(block, sz)
    @assert sz < block.sz "Cannot split a $block at too-large offset $sz"
    split = Block(block.ptr, sizeof(block) - sz; off = block.off + sz)
    block.sz = sz

    # update links
    split.prev = block
    split.next = block.next
    if block.next !== nothing
        block.next.prev = split
    end
    block.next = split

    return split
end

# merge a sequence of blocks that starts at `head` and ends with `tail` (inclusive)
function merge!(head, tail)
    head.next = tail.next
    if tail.next !== nothing
        tail.next.prev = head
    end
    head.sz = tail.sz + (tail.off - head.off)

    return head
end

@inline function actual_alloc(sz)
    ptr = CuArrays.actual_alloc(sz)
    block = ptr === nothing ? nothing : Block(ptr, sz)
end

function actual_free(block::Block)
    @assert iswhole(block) "Cannot free $block: block is not whole"
    if block.state != AVAILABLE
        error("Cannot free $block: block is not available")
    else
        CuArrays.actual_free(block.ptr)
        block.state = INVALID
    end
    return
end

## pooling

const pool_lock = ReentrantLock()

function scan!(pool, sz, max_overhead=typemax(Int))
    max_sz = max(sz + max_overhead, max_overhead)   # protect against overflow
    @lock pool_lock begin
        # get the first entry that is sufficiently large
        i = searchsortedfirst(pool, Block(nothing, sz; id=0))
        if i != pastendsemitoken(pool)
            block = deref((pool,i))
            @assert sizeof(block) >= sz
            if sz <= max_sz
                delete!((pool,i))   # FIXME: this allocates
                return block
            end
        end

        return nothing
    end
end

# compact sequences of blocks into larger ones.
# looks up possible sequences based on each block in the input set.
# destroys the input set.
# returns the net difference in amount of blocks.
function incremental_compact!(blocks)
    compacted = 0
    @lock pool_lock begin
        while !isempty(blocks)
            block = pop!(blocks)
            @assert block.state == AVAILABLE
            pool = get_pool(sizeof(block))
            @assert in(block, pool)

            # find the head of a sequence
            head = block
            while head.prev !== nothing && head.prev.state == AVAILABLE
                head = head.prev
                @assert in(head, pool)
            end
            szclass = size_class(sizeof(head))

            if head.next !== nothing && head.next.state == AVAILABLE
                delete!(pool, head)

                # find the tail (from the head, removing blocks as we go)
                tail = head.next
                while true
                    @assert szclass === size_class(sizeof(tail)) "block $tail should not have been split to a different pool than $head"
                    delete!(pool, tail)    # FIXME: allocates
                    delete!(blocks, tail)
                    tail.state = INVALID
                    compacted += 1
                    if tail.next !== nothing && tail.next.state == AVAILABLE
                        tail = tail.next
                    else
                        break
                    end
                end

                # compact
                head = merge!(head, tail)
                @assert !in(head, pool) "$head should not be in the pool"
                @assert head.state == AVAILABLE
                @assert szclass === size_class(sizeof(head)) "compacted $head should not end up in a different pool"
                push!(pool, head)
            end
        end
    end
    return compacted
end

function reclaim!(pool, sz=typemax(Int))
    freed = 0

    @lock pool_lock begin
        # mark non-split blocks
        candidates = Block[]
        for block in pool
            if iswhole(block)
                push!(candidates, block)
            end
        end

        # free them
        for block in candidates
            delete!(pool, block)
            freed += sizeof(block)
            actual_free(block)
            freed >= sz && break
        end
    end

    return freed
end

# repopulate the pools from the list of freed blocks
function repopulate()
    blocks = @safe_lock freed_lock begin
        isempty(freed) && return
        blocks = Set(freed)
        empty!(freed)
        blocks
    end

    @lock pool_lock begin
        for block in blocks
            pool = get_pool(sizeof(block))
            @assert !in(block, pool) "$block should not be in the pool"
            @assert block.state == FREED "$block should have been marked freed"
            block.state = AVAILABLE
            push!(pool, block) # FIXME: allocates
        end

        incremental_compact!(blocks)
    end

    return
end

const SMALL = 1
const LARGE = 2
const HUGE  = 3

# sorted containers need unique keys, which the size of a block isn't.
# mix in the block address to keep the key sortable, but unique.
unique_sizeof(block::Block) = (UInt128(sizeof(block))<<64) | UInt64(block.id)
const UniqueIncreasingSize = Base.By(unique_sizeof)

const pool_small = SortedSet{Block}(UniqueIncreasingSize)
const pool_large = SortedSet{Block}(UniqueIncreasingSize)
const pool_huge  = SortedSet{Block}(UniqueIncreasingSize)

const freed = Vector{Block}()
const freed_lock = SpinLock()

function size_class(sz)
    if sz <= SMALL_CUTOFF
        SMALL
    elseif SMALL_CUTOFF < sz <= LARGE_CUTOFF
        LARGE
    else
        HUGE
    end
end

@inline function get_pool(sz)
    szclass = size_class(sz)
    if szclass == SMALL
        return pool_small
    elseif szclass == LARGE
        return pool_large
    elseif szclass == HUGE
        return pool_huge
    end
end

function pool_alloc(sz)
    szclass = size_class(sz)

    # round off the block size
    req_sz = sz
    roundoff = if szclass == SMALL
        SMALL_ROUNDOFF
    elseif szclass == LARGE
        LARGE_ROUNDOFF
    elseif szclass == HUGE
        HUGE_ROUNDOFF
    end
    sz = cld(sz, roundoff) * roundoff
    szclass = size_class(sz)

    # select a pool
    pool = get_pool(sz)

    # determine the maximum scan overhead
    max_overhead = if szclass == SMALL
        SMALL_OVERHEAD
    elseif szclass == LARGE
        LARGE_OVERHEAD
    elseif szclass == HUGE
        HUGE_OVERHEAD
    end

    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 repopulate" repopulate()

        @pool_timeit "$phase.2 scan" begin
            block = scan!(pool, sz, max_overhead)
        end
        block === nothing || break

        @pool_timeit "$phase.3 alloc" begin
            block = actual_alloc(sz)
        end
        block === nothing || break

        # we're out of memory, try freeing up some memory. this is a fairly expensive
        # operation, so start with the largest pool that is likely to free up much memory
        # without requiring many calls to free.
        for pool in (pool_huge, pool_large, pool_small)
            @pool_timeit "$phase.4a reclaim" reclaim!(pool, sz)
            @pool_timeit "$phase.4b alloc" block = actual_alloc(sz)
            block === nothing || break
        end
        block === nothing || break

        # last-ditch effort, reclaim everything
        @pool_timeit "$phase.5a reclaim" begin
            reclaim!(pool_huge)
            reclaim!(pool_large)
            reclaim!(pool_small)
        end
        @pool_timeit "$phase.5b alloc" block = actual_alloc(sz)
    end

    if block !== nothing
        @assert size_class(sizeof(block)) == szclass "Cannot satisfy an allocation in pool $szclass with a buffer in pool $(size_class(sizeof(block)))"

        # split the block if that would yield one from the same pool
        # (i.e. don't split off chunks that would be way smaller)
        # TODO: avoid creating unaligned blocks here (doesn't currently seem to happen
        #       because of the roundoff, but we should take care anyway)
        remainder = sizeof(block) - sz
        if szclass != HUGE && remainder > 0 && size_class(remainder) == szclass
            split = split!(block, sz)
            split.state = AVAILABLE
            @lock pool_lock begin
                push!(pool, split)
            end
        end
    end

    return block
end

function pool_free(block)
    # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
    # and to simplify locking (preventing concurrent access during GC interventions)
    block.state = FREED
    @safe_lock freed_lock begin
        push!(freed, block)
    end
end

## interface

const allocated_lock = SpinLock()
const allocated = Dict{CuPtr{Nothing},Block}()

init() = return

function alloc(sz)
    block = pool_alloc(sz)
    if block !== nothing
        block.state = ALLOCATED
        ptr = pointer(block)
        @safe_lock allocated_lock begin
            @assert !haskey(allocated, ptr) "Newly-allocated block $block is already allocated"
            allocated[ptr] = block
        end
        return ptr
    else
        return nothing
    end
end

function free(ptr)
    block = @safe_lock allocated_lock begin
        block = allocated[ptr]
        delete!(allocated, ptr)
        block
    end
    block.state == ALLOCATED || error("Cannot free a $(block.state) block")
    pool_free(block)
    return
end

function reclaim(sz::Int=typemax(Int))
    repopulate()

    freed_sz = 0
    for pool in (pool_huge, pool_large, pool_small)
        freed_sz >= sz && break
        freed_sz += reclaim!(pool, sz-freed_sz)
    end
    return freed_sz
end

used_memory() = @safe_lock allocated_lock mapreduce(sizeof, +, values(allocated); init=0)

function cached_memory()
    sz = @safe_lock freed_lock mapreduce(sizeof, +, freed; init=0)
    @lock pool_lock for pool in (pool_small, pool_large, pool_huge)
        sz += mapreduce(sizeof, +, pool; init=0)
    end
    return sz
end

end


text:
using Test
include(joinpath(@__DIR__, "../src/ModelEnumeration.jl"))
include(joinpath(@__DIR__, "../data/sketches/cat.jl"))

F = catfls;

# Chase below explicitly  # chase_below(db, catfls, 1)
db = init_db(reset=false);

# There is one empty model
@test length(chase_below(db, catfls, 0)) == 1

# There is no model with only one O+A, since each O needs an identity A.
@test length(chase_below(db, catfls, 1)) == 1

# There are two models of size at most 2. The empty model, and the model with
# one (id) arrow and one object.
ms_2 = chase_below(db, catfls, 2)
@test length(ms_2) == 2
ms_2_again = chase_below(db, catfls, 2) # now results are already in DB
@test ms_2 == ms_2_again

# category with one object = monoid. There are two addn'l monoids with 2 arrows.
ms_3 = chase_below(db, catfls, 3)
@test length(ms_3) == 4

# Example from earlier
if 1+1==1
  include(joinpath(@__DIR__, "../data/sketches/fg.jl"))

  # chase a complete model, should give one model
  Ks, Ms = chasestep(fls, xmodrel)
  @test isempty(Ks)
  @test only(values(Ms)) == xmodrel

  db = init_db("test.db"; rem=true)
  add_model(db, fls, xmodrel) # adds one premodel and one model
  Kids, Mids = chasestep_db(db, fls, 1)
  @test isempty(Kids)
  @test Mids == [1]

  chase_below(db, fls, 1)
end

# I = grph_to_cset(:cat, catschema);
# add_part!(I, :A);
# add_part!(I, :O; refl=1);
# J = initrel(catfls, I);
# Ks = apply_tgds(catfls, J);
# K = Ks[1] # 3^3 models, first one is free
# e = apply_cones!(catfls, K)
# apply_egds!(catfls, K, e)
# Ks, Ms = chasestep(catfls, J)


text:
#=
This file defines additional functions to return objects of type Distribution. This is
necessary because we specify prior distributions wrt mean and SD
(for beta and gamma-distributed parameters) and ν and σ (for inverse gamma-distributed
parameters). Note these functions are NOT new methods for the Distributions.Beta, etc.
functions, but rather new functions with the same names.
=#

import Distributions: params, mean, std, pdf, logpdf, rand
import Base: length, rank

"""
```
BetaAlt(μ::AbstractFloat, σ::AbstractFloat)
```

Given μ and σ, calculate α and β and return a Distributions.Beta Distribution object.

### Arguments
`μ`: The mean of the desired distribution
`σ`: The standard deviation of the desired distribution
"""
function BetaAlt(μ::AbstractFloat, σ::AbstractFloat)
    α = (1-μ) * μ^2 / σ^2 - μ
    β = α * (1/μ - 1)
    return Distributions.Beta(α, β)
end

"""
```
GammaAlt(μ::AbstractFloat, σ::AbstractFloat)
```

Given μ and σ, calculate α and β and return a Distributions.Gamma object.

### Arguments
`μ`: The mean of the desired distribution
`σ`: The standard deviation of the desired distribution
"""
function GammaAlt(μ::AbstractFloat, σ::AbstractFloat)
    β = σ^2 / μ
    α = μ / β
    return Distributions.Gamma(α, β)
end

"""
```
type RootInverseGamma <: Distribution{Univariate, Continuous}
```

If x ~ RootInverseGamma(ν, τ²), then
  x² ~ ScaledInverseChiSquared(ν, τ²)
  x² ~ InverseGamma(ν/2, ντ²/2)
"""
type RootInverseGamma <: Distribution{Univariate, Continuous}
    ν::Float64
    τ::Float64
end

Distributions.params(d::RootInverseGamma) = (d.ν, d.τ)

"""
```
Distributions.pdf(d::RootInverseGamma, x::AbstractFloat)
```

Compute the pdf of a RootInverseGamma distribution at x.
"""
function Distributions.pdf(d::RootInverseGamma, x::AbstractFloat)
    (ν, τ) = params(d)
    return 2 * (ν*τ^2/2)^(ν/2) * exp((-ν*τ^2)/(2x^2)) / gamma(ν/2) / x^(ν+1)
end

"""
```
Distributions.logpdf(d::RootInverseGamma, x::AbstractFloat)
```

Compute the log pdf of a RootInverseGamma distribution at x.
"""
function Distributions.logpdf(d::RootInverseGamma, x::AbstractFloat)
    (ν, τ) = params(d)
    return log(2) - log(gamma(ν/2)) + (ν/2)*log(ν*τ^2/2) - ((ν+1)/2)*log(x^2) - ν*τ^2/(2x^2)
end

"""
```
DegenerateMvNormal <: Distribution{Multivariate, Continuous}
```

The DegenerateMvNormal type implements a degenerate multivariate normal distribution.
See [Multivariate normal distribution - Degenerate case](en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case).
"""
type DegenerateMvNormal <: Distribution{Multivariate, Continuous}
    μ::Vector          # mean
    σ::Matrix          # standard deviation
end

"""
```
rank(d::DegenerateMvNormal)
```

Returns the rank of d.σ.
"""
function rank(d::DegenerateMvNormal)
    return rank(d.σ)
end

"""
```
length(d::DegenerateMvNormal)
```

Returns the dimension of d.
"""
Base.length(d::DegenerateMvNormal) = length(d.μ)

"""
```
Distributions.rand{T<:AbstractFloat}(d::DegenerateMvNormal; cc::T = 1.0)
```

Generate a draw from d with variance optionally scaled by cc^2.
"""
function Distributions.rand{T<:AbstractFloat}(d::DegenerateMvNormal; cc::T = 1.0)
    return d.μ + cc*d.σ*randn(length(d))
end

"""
```
moments(dist::DegenerateMvNormal)
```

Compute the mean μ and standard deviation σ of a DSGE.DegenerateMvNormal object.
"""
function moments(dist::DegenerateMvNormal)
    return dist.μ, dist.σ
end

"""
```
moments(dist::RootInverseGamma)
```

Return the RootInverseGamma parameters τ and ν.
"""
function moments(dist::RootInverseGamma)
    ν, τ = params(dist)

    return τ, ν

    # Compute the mean μ and standard deviation σ of a DSGE.RootInverseGamma object.

    # A Root Inverse Gamma / Nagasaki Scaled Chi2 Distribution's mean and standard deviation
    # can be computed as follows:

    #     μ = √(β) * Γ(α - 0.5) / Γ(α)
    #     σ = √(β / (α - 1) - μ²)

    # where α = ν/2 and β = τ²*ν/2.

    # α = ν/2
    # β = τ^2*ν/2

    # μ = β^(.5) * gamma(α - 0.5) / gamma(α)
    # σ = (β / (α - 1) - μ^2)^(0.5)

    # return μ, σ
end

"""
```
moments(dist::Distributions.Beta)
```

Compute the mean μ and standard deviation σ of a Distributions.Beta object.
"""
function moments(dist::Distributions.Beta)
    α = dist.α
    β = dist.β

    μ = α / (α + β)
    σ = sqrt( α*β/((α+β)^2 * (α+β+1)) )
    return μ, σ
end

"""
```
moments(dist::Distributions.Gamma)
```

Compute the mean μ and standard deviation σ of a Distributions.Gamma object.
"""
function moments(dist::Distributions.Gamma)
    α = dist.α
    θ = dist.θ

    μ = α / θ
    σ = α / θ^2
    return μ, σ
end

"""
```
moments(dist::Distributions.Normal)
```

Return the mean and std dev or a Distributions.Normal object.
"""
function moments(dist::Distributions.Normal)
    params(dist)
end


text:
module PlottingModule

using Plots
using Printf

using ..ObjectAgentModule
using ..ActuatorArrayModule

export plot_AA, visu_pressure

# Plotting function
@recipe function plot_AA(aa::ActuatorArray{T,U};
    actuatorCommands = fill(NaN, aa.nx, aa.ny),
    actuatorType = :currents, #:phases
    boxes = Array{Array{Tuple{U,U},1},1}(),
    agents = Array{ObjectAgent, 1}(),
    agents_clrs = missing,
    dev_forces = Array{Tuple{T,T,T}, 1}(),
    aux_forces = Array{Tuple{T,T,T}, 1}(),
    ref_pos = Array{Tuple{T,T,T}, 1}(),
    distCircles = missing,
    Fscale = 1,
    plotUsedActuators = true,
    plotUsedActuators_enable = missing,
    showAgentNames = true,
    showReqForces = false,
    lwidth = 8,
    Fdev_clr = :yellow,
    Faux_clr = :green,
    Freq_clr = :red)  where {T<:Real,U<:Unsigned}
    # if actuator commands are provided, check whether it has the correct size
    @assert size(actuatorCommands) == (aa.nx, aa.ny) "The matrix of phase-shifts doesn't have the correct size!"

    if plotUsedActuators && length(agents) > 0 && ismissing(plotUsedActuators_enable)
        plotUsedActuators_enable = fill(true, length(agents));
    end

    aspect_ratio := :equal
    framestyle := :border

    currents_colorfun(α) = begin
        RGBA(1, 1-α, 1-α, 1);
    end
    phases_colorfun(α) = begin
        α = mod(α, 2*π)/2/π

        if α < 1/3
            k = 3*α
            clr = RGBA(1,1-k,1-k,1)
        elseif α < 2/3
            k = 3*(α-1/3)
            clr = RGBA(1-k,0,k,1)
        else
            k = 3*(α-2/3)
            clr = RGBA(k,k,1,1)
        end

        clr
    end

    if actuatorType == :currents
        act_colorfun = currents_colorfun
    elseif actuatorType == :phases
        act_colorfun = phases_colorfun
    else
        error("Unsupported actuator type.")
    end

    xlims := (-aa.dx/2, (aa.nx-0.5)*aa.dx)
    ylims := (-aa.dx/2, (aa.ny-0.5)*aa.dx)

    ticks := nothing

    # Generate boxes for plotting the used actuators
    x_box = [-aa.dx/2, aa.dx/2, aa.dx/2, -aa.dx/2];
    y_box = [-aa.dx/2, -aa.dx/2, aa.dx/2, aa.dx/2];

    # if boxes are provided plot the boxes around the actuators
    if length(boxes) != 0
        box_clrs = distinguishable_colors(length(boxes), lchoices = range(40, stop=80, length=15), cchoices = range(30, stop=100, length=15))
        k = 1;
        for k in eachindex(boxes)
            for act_box in boxes[k]
                @series begin
                    xc, yc = actuatorPosition(aa, act_box)

                    fill  := (0, box_clrs[k]);
                    label := "";
                    linewidth := 0;

                    x_box .+ xc, y_box .+ yc
                end
            end
        end
    end

    if length(agents) != 0
        if ismissing(agents_clrs)
            agents_clrs = distinguishable_colors(length(agents), lchoices = range(40, stop=80, length=15), cchoices = range(30, stop=80, length=15));
        end
    end

    # If a list of agents is provided, plot boxes around actuators used by the agents
    if length(agents) != 0 && plotUsedActuators
        for (k, agent) in enumerate(agents)
            if !plotUsedActuators_enable[k]
                continue;
            end

            if ismissing(distCircles)
                for (ak, act_box) in enumerate(agent.actList)
                    if ak ∉ agent.act_used
                        clralpha = 0.25;
                    else
                        clralpha = 0.5;
                    end
                    @series begin
                        xc, yc = actuatorPosition(aa, act_box)

                        fill  := (0, clralpha, agents_clrs[k])
                        label := ""
                        linealpha := 0

                        x_box .+ xc, y_box .+ yc
                    end
                end
            else
                th = range(0, 2*π, length=25);
                xc, yc = agent.pos;
                @series begin
                    fill  := (0, 0.5, agents_clrs[k])
                    seriescolor := agents_clrs[k]
                    label := ""
                    linealpha := 1

                    sin.(th).*distCircles[1] .+ xc, cos.(th).*distCircles[1] .+ yc
                end
                @series begin
                    fill  := (0, 0.25, agents_clrs[k])
                    seriescolor := agents_clrs[k]
                    label := ""
                    linealpha := 1

                    sin.(th).*distCircles[2] .+ xc, cos.(th).*distCircles[2] .+ yc
                end
            end
        end
    end

    # Generate circles for plotting the actuators
    if aa.shape == :circle || aa.shape == :coil
        th = range(0, 2*π, length=25);
        x_act_box = aa.width*cos.(th)./2;
        y_act_box = aa.width*sin.(th)./2;
    elseif aa.shape == :square
        x_act_box = [-aa.width/2, aa.width/2, aa.width/2, -aa.width/2, -aa.width/2];
        y_act_box = [-aa.width/2, -aa.width/2, aa.width/2, aa.width/2, -aa.width/2];
    else
        error("Unsupported shape of the actuators. The shape must be either :circle or :square")
    end

    # Iterate over indexes of the actuators
    for (ai, aj) in aa
        # skip the actuator if it is in the ignore list
        if in((ai, aj), aa.ignore_list)
            continue;
        end

        # get the position of the acutator with indexes (ai, aj)
        ax, ay = actuatorPosition(aa, (ai, aj));

        @series begin
            seriescolor := :black
            label := ""

            if ~isnan(actuatorCommands[ai,aj])
                fill_clr = act_colorfun(actuatorCommands[ai,aj])
            else
                fill_clr = RGBA(1, 1, 1, .8);
            end
            # fill  := (0, fill_clr)
            fill  := (0, fill_clr)

            x_act_box .+ ax, y_act_box .+ ay
        end

        if aa.shape == :coil
            @series begin
                seriescolor := :black
                label := ""

                fill  := (0, RGBA(0.6, 0.6, 0.6, 1))

                0.5.*x_act_box .+ ax, 0.5.*y_act_box .+ ay
            end
        end
    end

    # If a list of agents is provided, plot dots at the places of the agents
    if length(agents) != 0
        if !isa(agents[1], ObjectAgent_MAG) && !isa(agents[1], ObjectAgent_DEP)
            for (k, agent) in enumerate(agents)
                @series begin
                    markershape := :circle;
                    markersize := 10;
                    markercolor := agents_clrs[k];
                    if showAgentNames
                        label := agent.name;
                    else
                        label := "";
                    end
                    legendfontsize := 4;
                    linealpha := 0;

                    xc, yc = agent.pos;
                    [xc], [yc]
                end
            end
        else
            if isa(agents[1], ObjectAgent_MAG)
                rad = 10e-3;
            else
                rad = 25e-6;
            end

            th = range(0, 2*π, length=25);
            x_agt_box = rad.*cos.(th);
            y_agt_box = rad.*sin.(th);

            for (k, agent) in enumerate(agents)
                xc, yc = agent.pos;
                @series begin
                    seriescolor := :black
                    fill  := (0, agents_clrs[k]);
                    if showAgentNames
                        label := agent.name;
                    else
                        label := "";
                    end
                    legendfontsize := 4;

                    xc .+ x_agt_box, yc .+ y_agt_box
                end
            end
        end

    end

    # If agents and forces are provided, visualize force at the locations of the agents (Works only for ObjectAgent_DEP)
    if length(agents) != 0 && length(dev_forces) != 0 || showReqForces
        for (k, agent) in enumerate(agents)
            @assert (agent isa ObjectAgent_DEP) || (agent isa ObjectAgent_MAG) "Forces can be plotted only for ObjectAgent_DEP or ObjectAgent_MAG"

            xc, yc = agent.pos;
            # Developed force (provided in dev_forces[])
            if length(dev_forces) != 0
                @series begin
                    linecolor := Fdev_clr;
                    # linestyle := :dot;
                    linewidth := lwidth;
                    label := "";

                    [xc, xc + Fscale*dev_forces[k][1]*aa.dx], [yc, yc + Fscale*dev_forces[k][2]*aa.dx]
                end
            end

            # If available, plot also forces in  aux_forces[]
            if length(aux_forces) != 0
                @series begin
                    linecolor := Faux_clr;
                    # linestyle := :dot;
                    linewidth := lwidth;
                    label := "";

                    [xc, xc + Fscale*aux_forces[k][1]*aa.dx], [yc, yc + Fscale*aux_forces[k][2]*aa.dx]
                end
            end

            # Required force (stored by each agents in agent.Fdes)
            if showReqForces
                @series begin
                    linecolor := Freq_clr;
                    # linestyle := :dash;
                    label := "";
                    linealpha := 0.65;
                    linewidth := 0.65*lwidth;

                    [xc, xc + Fscale*agent.Fdes[1]*agent.Fdes_sc*aa.dx], [yc, yc + Fscale*agent.Fdes[2]*agent.Fdes_sc*aa.dx]
                end
            end

            if length(agent.Fdes) == 3 && showReqForces && length(dev_forces) != 0
                # Visualize relativr error in z-component of the developed force
                Fz_relerr = abs((agent.Fdes[3]*agent.Fdes_sc - dev_forces[k][3])/agent.Fdes[3]);
                @series begin
                    markershape := :circle;
                    markersize := 10*Fz_relerr;
                    markercolor := :red;
                    label := "";
                    linewidth := 0;
                    markerstrokewidth := 0;

                    [xc], [yc]
                end
            end
        end
    end

    # If reference positions are provided, plot these
    if length(agents) != 0 && length(ref_pos) != 0
        for (k, agent) in enumerate(agents)
            @series begin
                x_ref, y_ref = ref_pos[k];
                markershape := :circle;
                markersize := 10;
                markercolor := agents_clrs[k];
                markeralpha := 0.75;
                markerstrokewidth := 4;
                markerstrokealpha := 0.75;
                markerstrokecolor := :white;
                linewidth := 0;
                label := "";

                [x_ref], [y_ref]
            end
        end
    end
end

@recipe function visu_pressure(aa::ActuatorArray{T,U},
    z0::T,
    phases::Array{T,2};
    agents = Array{ObjectAgent, 1}(),
    N = 75) where {T<:Real, U<:Unsigned}
    @assert size(phases) == (aa.nx, aa.ny) "The matrix of phase-shifts doesn't have the correct size!"

    aspect_ratio := :equal
    framestyle := :border
    xlims := (aa.dx/2, (aa.nx-0.5)*aa.dx)
    ylims := (aa.dx/2, (aa.ny-0.5)*aa.dx)

    ticks := nothing

    xv = LinRange(-aa.dx/2,aa.nx*aa.dx-aa.dx/2, N);
    yv = LinRange(-aa.dx/2,aa.nx*aa.dx-aa.dx/2, N);
    P = pressureField(aa, z0, phases, xv, yv);

    @series begin
        seriestype := :heatmap;
        xv, yv, P
    end

    # If a list of agents is provided, plot dots at the places of the agents
    if length(agents) != 0
        agents_clrs = distinguishable_colors(length(agents), lchoices = range(40, stop=80, length=15), cchoices = range(30, stop=80, length=15))
        for (k, agent) in enumerate(agents)
            @series begin
                markershape := :circle;
                markersize := 4;
                markercolor := agents_clrs[k];
                # label := agent.name;
                # legendfontsize := 4;
                label := "";
                linealpha := 0;
                seriestype := :path;

                markerstrokecolor := :black;
                markerstrokewidth := 1;

                xc, yc, _ = agent.pos;
                [xc], [yc]
            end
        end
    end
end

end


text:
@testset "ECDSA" begin

    import Secp256k1: Point, G, N

    @testset "Signature Verification" begin
        kp = KeyPair{:ECDSA}(rand(big.(1:N-1)))
        𝑧  = rand(Int)
        𝑠  = ECDSA.sign(kp, 𝑧)
        @test ECDSA.verify(kp.𝑄, 𝑧, 𝑠)

        point = Point(
            big"0x887387e452b8eacc4acfde10d9aaf7f6d9a0f975aabb10d006e4da568744d06c",
            big"0x61de6d95231cd89026e286df3b6ae4a894a3378e393e93a0f45b666329a0ae34")
        z = big"0xec208baa0fc1c19f708a9ca96fdeff3ac3f230bb4a7ba4aede4942ad003c0f60"
        r = big"0xac8d1c87e51d0d441be8b3dd5b05c8795b48875dffe00b7ffcfac23010d3a395"
        s = big"0x68342ceff8935ededd102dd876ffd6ba72d6a427a3edb13d26eb0781cb423c4"
        @test ECDSA.verify(point, z, Signature{:ECDSA}(r, s))

        z = big"0x7c076ff316692a3d7eb3c3bb0f8b1488cf72e1afcd929e29307032997a838a3d"
        r = big"0xeff69ef2b1bd93a66ed5219add4fb51e11a840f404876325a1e8ffe0529a2c"
        s = big"0xc7207fee197d27c618aea621406f6bf5ef6fca38681d82b2f06fddbdce6feab6"
        @test ECDSA.verify(point, z, Signature{:ECDSA}(r, s))
    end
end


text:
using Plots

# n=100; xr= 10*rand(n); yr= 10*rand(n); zr= Float64[ sin(xr[i])+cos(yr[i]) for i=1:n];
# surface(xr, yr, zr, size=[800,480])

# x = repeat(0:2/20:2,21)
# y = reshape(repeat(0:2/20:2,21),21,21)'[:]
# f(x,y) = @. sin(x)+cos(y)
# surface(x,y,f(x,y),colorbar=false, fc = :haline, legend=false)

x= [1.0:0.1:5;]; y= [4.0:0.1:9.0;]; f(x,y)= sin(x)*cos(y)+4;
p1= surface(x,y,f,camera=(30,30),legend=false,size=[800,500],xlims=[-1,6],zlims=[2,5],fc =:haline,grid=true,axis=false,grid_linewidth=3);	# 30 azimuth, 30 elevation is the default
plot(p1)

savefig("3d-surface.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
#=
You are given an M by N matrix consisting of booleans that represents a board. Each True boolean represents a wall. Each False boolean represents a tile you can walk on.

Given this matrix, a start coordinate, and an end coordinate, return the minimum number of steps required to reach the end coordinate from the start. If there is no possible path, then return null. You can move up, left, down, and right. You cannot move through walls. You cannot wrap around the edges of the board.

For example, given the following board:

[[f, f, f, f],
[t, t, f, t],
[f, f, f, f],
[f, f, f, f]]

and start = (3, 0) (bottom left) and end = (0, 0) (top left), the minimum number of steps required to reach the end is 7, since we would need to go through (1, 2) because there is a wall everywhere else on the second row.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using CSVFiles
using DataFrames
using DrWatson

ddir = projectdir("figures", "figure_gene_expr", "data")

comb_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-10-04.csv")))
dmc_data = comb_data[comb_data.Alg.=="Direct MC", :]
smc_data = comb_data[comb_data.Alg.=="SMC", :]
smc_estimate = groupby(smc_data, :M)[(; M = 128)]
ti_data = comb_data[(comb_data.Alg.=="TI").&(comb_data.Duration.==comb_data.DiscreteTimes), :]
zechner_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-18_zechner.csv")))
gaussian_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-07-19_gaussian.csv")))

using CairoMakie

# presentation_theme = Theme(font = "Noto Sans Regular", linewidth=3)

# fig1 = with_theme(presentation_theme) do
#     f = Figure(resolution=(600, 800))
#     ax = Axis(f[1, 1], xlabel="trajectory duration", ylabel="mutual information (nats)")
#     dmc_d = groupby(dmc_data, :M)[end]
#     dmc_lines = lines!(ax, dmc_d.DiscreteTimes, dmc_d.Mean, label="direct")
#     smc_lines = lines!(ax, smc_estimate.DiscreteTimes, smc_estimate.Mean, color="black")
#     zec_lines = lines!(zechner_data.Duration, zechner_data.PMI, linewidth=2, linestyle=:dash, color=:red)

#     legend = Legend(
#         f[2, 1],
#         [dmc_lines, smc_lines, zec_lines],
#         ["direct", "segment-by-segment", "Duso, Zechner (2019)"],
#         orientation = :horizontal,
#         tellheight = true
#         # [group_color, ti_color, smc_color, approx_entry],
#         # [m_label, ["256 MCMC steps", "4096 MCMC steps"],
#         # ["M = 128"], ["Duso, Zechner (2019)", "Gaussian"]],
#         # ["brute force", "thermodynamic\nintegration", "particle filter", "approximations"],
#         # rowgap=0,
#         # titlegap=2,
#         # groupgap=5,
#         # patchsize=(10,10),
#         # padding=(3,3,3,3),
#         # linewidth=1.0,
#         # gridshalign=:left,
#         # titlehalign=:left,
#         # titlefont="Noto Sans Medium",
#         # framevisible=false
#     )
#     xlims!(ax, 0, 8)
#     ylims!(ax, 0, 20)

#     f
# end
# save(projectdir("plots", "figure_gene_expr", "gene-expr-figure-p.pdf"), fig1, pt_per_unit = 1)

normalfont = projectdir("figures", "fonts", "NotoSans-Condensed.ttf")
titlefont = projectdir("figures", "fonts", "NotoSans-SemiBoldItalic.ttf")

fontsize_theme = Theme(
    font = normalfont,
    fontsize = 8,
    colgap = 3,
    linewidth = 2,
    markersize = 6,
    spinewidth = 0.5,
    Axis = (
        spinewidth = 1.0,
        xgridwidth = 1.0,
        ygridwidth = 1.0,
        xtickwidth = 1.0,
        ytickwidth = 1.0,
        xticksize = 3,
        yticksize = 3,
        xlabelpadding = 0,
        ylabelpadding = 5,
        xticklabelpad = 0,
        yticklabelpad = 3,
    ),
    Legend = (
        framewidth = 1.0,
        rowgap = 0,
        titlegap = 2,
        groupgap = 5,
        patchsize = (11, 10),
        padding = (3, 3, 3, 3),
        linewidth = 1.25,
        gridshalign = :left,
        titlehalign = :left,
        titlefont = titlefont,
        framevisible = false
    )
)

fig = with_theme(fontsize_theme) do
    f = Figure(resolution = (246, 1.4 * 146), figure_padding = 2)
    Axis(f[1, 1], xlabel = "trajectory duration", ylabel = "mutual information (nats)")

    groups = collect(groupby(dmc_data, :M; sort = true))

    m_label = ["M = $(g.M[1])" for g in groups]
    cmap = cgrad(:linear_kgy_5_95_c69_n256, rev = true)
    colors = [cmap[i/length(groups)*0.7] for i = 1:length(groups)]
    ti_colors = [cgrad(:lajolla)[i] for i = range(0.2, 0.6, length = 3)]
    pathmi_color = :orangered

    for (i, group) in enumerate(groups)
        label = "$(group.M[1])"
        lines!(group.DiscreteTimes, group.Mean, linewidth = 0.75, label = label, color = colors[i])
    end

    group_color = [LineElement(color = color) for color in colors]
    smc_color = [LineElement(color = :black)]

    lines!(smc_estimate.DiscreteTimes, smc_estimate.Mean, label = "SMC", linewidth = 1.7, color = :black)
    lines!(zechner_data.Duration, zechner_data.PMI, linewidth = 1.5, linestyle = :dash, color = pathmi_color)
    lines!(gaussian_data.DiscreteTimes, gaussian_data.Value, linewidth = 1.5, linestyle = :dot, color = "#2b8cbe")

    ti_label = []
    ti_color = []
    for (i, group) in enumerate(groupby(ti_data, :M, sort = true))
        scatter!(group.DiscreteTimes, group.Mean, marker = :cross, color = ti_colors[i])
        push!(ti_label, "M = $(group.M[1])")
        push!(ti_color, MarkerElement(marker = :cross, color = ti_colors[i]))
    end

    approx_entry = [LineElement(linestyle = :dash, color = pathmi_color), LineElement(linestyle = :dot, color = "#2b8cbe")]

    legend = Legend(
        f[1, 2],
        [group_color, ti_color, smc_color, approx_entry],
        [m_label, ti_label,
            ["M = 128"], ["Duso, et. al. (2019)", "Gaussian (LNA)"]],
        ["DPWS", "TI-PWS", "RR-PWS", "approximations"],
    )
    xlims!(0.0, 10.2)
    ylims!(0.0, 42)

    inset_lims = (3.3, 8.5)
    lines!([(0, 0), (inset_lims[1], 0), inset_lims, (0, inset_lims[2]), (0, 0)], color = :black, linewidth = 0.75, linestyle = :dot)

    inset_ax = Axis(f[1, 1],
        width = Relative(0.55),
        height = Relative(0.33),
        halign = 0.25,
        valign = 0.95,
        xtickalign = 1.0,
        ytickalign = 1.0,
        backgroundcolor = :white,
        xticksize = 2,
        xticklabelsize = 7,
        xgridvisible = true,
        yticksize = 2,
        yticklabelsize = 7,
    )

    for (i, group) in enumerate(groups)
        label = "$(group.M[1])"
        lines!(inset_ax, group.DiscreteTimes, group.Mean, linewidth = 0.5, label = label, color = colors[i])
    end
    lines!(inset_ax, smc_estimate.DiscreteTimes, smc_estimate.Mean, linewidth = 1.0, color = :black)
    lines!(inset_ax, zechner_data.Duration, zechner_data.PMI, linewidth = 0.75, linestyle = :dash, color = pathmi_color)
    lines!(inset_ax, gaussian_data.DiscreteTimes, gaussian_data.Value, linewidth = 1.5, linestyle = :dot, color = "#2b8cbe")

    for (i, group) in enumerate(groupby(ti_data, :M, sort = true))
        scatter!(inset_ax, group.DiscreteTimes, group.Mean, marker = :cross, color = ti_colors[i])
    end

    xlims!(inset_ax, 0, inset_lims[1])
    ylims!(inset_ax, 0, inset_lims[2])
    translate!(inset_ax.scene, 0, 0, 10)
    translate!(inset_ax.elements[:background], 0, 0, 9)
    translate!(inset_ax.elements[:xgridlines], 0, 0, 9)
    translate!(inset_ax.elements[:ygridlines], 0, 0, 9)

    f
end

save(projectdir("figures", "figure_gene_expr", "gene-expr-figure.svg"), fig, pt_per_unit = 1)
save(projectdir("figures", "figure_gene_expr", "gene-expr-figure.pdf"), fig, pt_per_unit = 1)

fig