text:
# Semidiscrete KT Scheme: Second-Order
# Based on:
# Kurganov, Tadmor. New High Resolution Central Schemes for Non Linear Conservation
# Laws and Convection-Difussion Equations. Journal of Comp Physics 160, pp 241-282. 2000

struct FVSKTScheme{ltype <: AbstractSlopeLimiter} <: AbstractFVAlgorithm
  slopeLimiter :: ltype
end

function FVSKTScheme(;slopeLimiter=GeneralizedMinmodLimiter())
  FVSKTScheme(slopeLimiter)
end

function update_flux_value(u,∇u,j,mesh,Flux, alg::FVSKTScheme)
    # Local speeds of propagation
    uminus=cellval_at_left(j,u,mesh)+0.5*cellval_at_left(j,∇u,mesh)
    uplus=cellval_at_right(j,u,mesh)-0.5*cellval_at_right(j,∇u,mesh)
    ul = cellval_at_left(j,u,mesh)
    ur = cellval_at_right(j,u,mesh)
    aa = max(fluxρ(uminus,Flux),fluxρ(uplus,Flux))
    # Numerical Fluxes
    0.5*(Flux(uplus)+Flux(uminus)) - aa/2*(uplus - uminus)
end

"""
compute_fluxes!(hh, Flux, u, mesh, dt, alg::FVSKTScheme, ::Type{Val{true}})
Numerical flux of Kurkanov Tadmor scheme in 1D
"""
function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::FVSKTScheme, noscalar::Bool, ::Type{Val{true}})
    slopeLimiter = alg.slopeLimiter
    #update vector
    # 1. slopes
    ∇u = compute_slopes(u, mesh, slopeLimiter, noscalar, Val{true})

    #update vector
    for j in node_indices(mesh)
        if noscalar
            fluxes[:,j] .= update_flux_value(u,∇u,j,mesh,Flux,alg)
        else
            fluxes[j] = update_flux_value(u,∇u,j,mesh,Flux,alg)
        end
    end
end

function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::FVSKTScheme, noscalar::Bool, ::Type{Val{false}})
    slopeLimiter = alg.slopeLimiter
    #update vector
    # 1. slopes
    ∇u = compute_slopes(u, mesh, slopeLimiter, noscalar, Val{false})

    #update vector
    for j in node_indices(mesh)
        if noscalar
            fluxes[:,j] .= update_flux_value(u,∇u,j,mesh,Flux,alg)
        else
            fluxes[j] = update_flux_value(u,∇u,j,mesh,Flux,alg)
        end
    end
end


text:
"""
    InfExtendedReal{T} <: Real

The type `T` extended with positive and negative infinity.
"""
struct InfExtendedReal{T<:Real} <: Real
  flag :: InfFlag
  finitevalue :: T

  InfExtendedReal{T}(x::T) where {T<:Real} = new(FINITE, x)
  InfExtendedReal{T}(x::Infinite) where {T<:Real} = new(x==PosInf ? POSINF : NEGINF)

  # For internal testing use only
  function Base.reinterpret(::Type{InfExtendedReal{T}}, t::Tuple{InfFlag, T}) where T<:Real
      new{T}(t[1], t[2])
  end
end

# Since InfExtendedReal is a subtype of Real, and Infinite is also a subtype of real,
# we can just use `x.val` to get either the finite value, or the infinite value. This will
# make arithmetic much simpler.
function Base.getproperty(x::InfExtendedReal, s::Symbol)
    if s === :val
        return x.flag != FINITE ? (x.flag == NEGINF ? -∞ : ∞) : x.finitevalue
    else
        return getfield(x, s)
    end
end

InfExtendedReal{T}(x::Real) where {T<:Real} = InfExtendedReal{T}(isinf(x) ? convert(Infinite, x) : convert(T, x))
InfExtendedReal{T}(x::InfExtendedReal) where {T<:Real} = InfExtendedReal{T}(x.val)

"""
    InfExtendedReal(T)

The union of `T` and `Infinite`: either `T` if infinity can be represented as a `T`, or else `InfExtendedReal{T}`.
"""
@generated InfExtendedReal(::Type{T}) where {T<:Real} = hasinf(T) ? T : InfExtendedReal{T}

"""
    InfExtendedReal(x)

Converts `x` to a `InfExtendedReal(typeof(x))`.
"""
@generated InfExtendedReal(x::T) where {T<:Real} = hasinf(T) ? :x : :($(InfExtendedReal(T))(x))

Utils.isposinf(x::InfExtendedReal) = x.flag == POSINF || x.flag == FINITE && isposinf(x.finitevalue)
Utils.isneginf(x::InfExtendedReal) = x.flag == NEGINF || x.flag == FINITE && isneginf(x.finitevalue)

ext: jl


text:
using JuMP
using Combinatorics

# Inspired from Joey Huchette's test in ConvexHull.jl
function hypercubetest(lib::PolyhedraLibrary, n)
    m = Model()
    @variable(m, 0 ≤ x[1:n] ≤ 1)

    poly = polyhedron(m, lib)

    row = 0
    V = zeros(Int, 2^n, n)
    for k in 0:n
        for p in combinations(1:n, k)
            row += 1
            V[row, p] = ones(Int, length(p))
        end
    end
    generator_fulltest(poly, V)
end


text:
export chebyshev, chebyshev!

chebyshev(A, b, λmin::Real, λmax::Real, Pr = 1, n = size(A,2);
          tol::Real = sqrt(eps(typeof(real(b[1])))), maxiter::Int = n^3) =
    chebyshev!(zerox(A, b), A, b, λmin, λmax, Pr, n; tol=tol,maxiter=maxiter)

function chebyshev!(x, A, b, λmin::Real, λmax::Real, Pr = 1, n = size(A,2);
	tol::Real = sqrt(eps(typeof(real(b[1])))), maxiter::Int = n^3
    )
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm,maxiter)

	K = KrylovSubspace(A, n, 1, Adivtype(A, b))
	init!(K, x)
	chebyshev_method!(history, x, K, b, λmin, λmax, Pr; tol = tol, maxiter = maxiter)
    x, history
end

function chebyshev_method!(
    log::ConvergenceHistory, x, K::KrylovSubspace, b, λmin::Real, λmax::Real,
    Pr = 1; tol::Real = sqrt(eps(typeof(real(b[1])))), maxiter::Int = K.n^3
    )

	local α, p
	K.order = 1
    tol = tol*norm(b)
	r = b - nextvec(K)
	d::eltype(b) = (λmax + λmin)/2
	c::eltype(b) = (λmax - λmin)/2
	for iter = 1:maxiter
        nextiter!(log, mvps=1)
		z = Pr\r
		if iter == 1
			p = z
			α = 2/d
		else
			β = (c*α/2)^2
			α = 1/(d - β)
			p = z + β*p
		end
		append!(K, p)
		update!(x, α, p)
		r -= α*nextvec(K)
		#Check convergence
		resnorm = norm(r)
        push!(log, :resnorm, resnorm)
        resnorm < tol && break
	end
    shrink!(log)
    setconv(log, 0<=norm(r)<tol)
	x
end


text:
module VarianceComponentsHDFE

# using DataDeps
using CSV
using DataFrames, Parameters
using LinearAlgebra, SparseArrays, Random, Statistics
using SparseArrays, LightGraphs
using Distributions
using LinearOperators, FastClosures, Krylov
using ArgParse
using AlgebraicMultigrid

#include("init.jl")
include("leave_out_correction.jl")
include("parameters_settings.jl")
include("laplacians/Laplacians.jl")

using .Laplacians
include("solvers.jl")

export find_connected_set,prunning_connected_set,drop_single_obs
export compute_movers, accumarray
export lincom_KSS, compute_matchid, leave_out_estimation, get_leave_one_out_set
export VCHDFESettings, JLAAlgorithm, ExactAlgorithm, AbstractLeverageAlgorithm
export leave_out_KSS, leverages, kss_quadratic_form, sigma_for_stayers

# Exporting these functions for ease of benchmarking/testing
export computeLDLinv, approxcholOperator, approxcholSolver

function parse_commandline()
    argparsesettings_obj = ArgParseSettings()

    @add_arg_table! argparsesettings_obj begin
        "path"
            help = "path to CSV file containing data"
            required = true
        "--first_id"
            help = "column index in CSV file for the first ID (e.g. Person).  Use the most granular type."
            arg_type = Int
            default = 1
        "--second_id"
            help = "column index in CSV file for the second ID (e.g. Firm).  Use the less granular type."
            arg_type = Int
            default = 2
        "--outcome_id"
            help = "column index in CSV file for outcome (e.g. Wage)."
            arg_type = Int
            default = 4
        "--no_first_id_effects"
            help = "No computing and showing of first_id effects"
            action = :store_true
        "--no_cov_effects"
            help = "No computing and showing of covariace effects"
            action = :store_true
        "--leave_out_level"
            help = "leave out level: obs or match. It defaults to match."
            arg_type = String
            default = "match"
        "--algorithm"
            help = "type of algorithm: Exact or JLA. It defaults to be Exact if the number of observations is less than 5000, and JLA otherwise."
            arg_type = String
            default = "Default"
        "--simulations"
            help = "number of simulations in the JLA algorithm. If 0, defaults to 200"
            arg_type = Int
            default = 0
        "--header"
            help = "CSV file contains header"
            action = :store_true
        "--first_id_display"
            help = "The display text associated with first_id (e.g. Person)."
            arg_type = String
            default = "Person"
        "--second_id_display"
            help = "The display text associated with second_id (e.g. Firm)"
            arg_type = String
            default = "Firm"
        "--outcome_id_display"
            help = "The display text associated with outcome_id (e.g. Wage)"
            arg_type = String
            default = "Wage"
        "--detailed_csv_path"
            help = "path to the CSV for the detailed output for each observable"
            arg_type = String
            default = joinpath(pwd(), "variance_components.csv")
        "--results_path"
            help = "path to the results of the output"
            arg_type = String
            default = joinpath(pwd(), "results.txt")
        "--write_detailed_csv"
            help = "write the detailed output to a CSV"
            action = :store_true
        "--write_results"
            help = "write the results to a file"
            action = :store_true
        #adding the print level command
        "--print_level"
            help = "Level of verbosity of output."
            arg_type = Int
            default = 1
        "--do_lincom"
            help = "perform linear combination based on observables Z"
            action = :store_true
        #Select the number of columns or put labels
        "--covariates"
            help = "Covariates. You can set column numbers or String name."
            nargs = '*'
            default = []
        #Select the number of columns or put labels
        "--lincom_covariates"
            help = "Lincom covariates. You can set column numbers or String name."
            nargs = '*'
            default = []
    end

    return parse_args(argparsesettings_obj)
end

function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    parsed_args = parse_commandline()

    println("Present working directory is $(pwd())")

    path = parsed_args["path"]
    #TODO: not sure about the paths
    if isabspath(path) == false
        path = joinpath(pwd(), path)
    end
    header = parsed_args["header"]
    first_idx = parsed_args["first_id"]
    second_idx = parsed_args["second_id"]
    outcome_idx = parsed_args["outcome_id"]
    leave_out_level = parsed_args["leave_out_level"]
    algorithm = parsed_args["algorithm"]
    first_id_effects = parsed_args["no_first_id_effects"] == 1 ? 0 : 1
    cov_effects = parsed_args["no_cov_effects"] == 1 ? 0 : 1
    simulations = parsed_args["simulations"]
    first_id_display = parsed_args["first_id_display"]
    second_id_display = parsed_args["second_id_display"]
    outcome_id_display = parsed_args["outcome_id_display"]
    print_level = parsed_args["print_level"]
    covariates = parsed_args["covariates"]
    lincom_covariates = parsed_args["lincom_covariates"]

    algorithm = lowercase(algorithm)
    leave_out_level = lowercase(leave_out_level)
    first_id_display_small = lowercase(first_id_display)
    second_id_display_small = lowercase(second_id_display)

    print_level > 0 && println("Number of threads: $(Threads.nthreads())")

    data  = DataFrame(CSV.File(path; header=header))
    first_id = data[:,first_idx]
    second_id = data[:,second_idx]
    y = data[:,outcome_idx]

    if covariates == []
        controls = nothing
    else
        #It creates a dictionary of the covariates
        if typeof(covariates[1]) == String
            #Build controls matrix
            controls = data[:,covariates[1]]
            if length(covariates)>=2
                for k=2:length(covariates)
                    hcat(controls, data[:,covariates[k]])
                end
            end
        else
            println("WARNING: Elements of covariates are not defined correctly. No partialling out will be performed.")
            controls = nothing
        end
    end

    if algorithm == "default"
        if length(y) > 5000
            algorithm = "jla"
        else
            algorithm = "exact"
        end
    end

    if algorithm == "exact"
        #todo work with settings arguments
        settings = VCHDFESettings(leverage_algorithm =  ExactAlgorithm(),
            first_id_effects =first_id_effects,
            cov_effects = cov_effects,
            first_id_display = first_id_display,
            first_id_display_small = lowercase(first_id_display),
            second_id_display = second_id_display,
            second_id_display_small = lowercase(second_id_display),
            outcome_id_display= outcome_id_display,
            outcome_id_display_small = lowercase(outcome_id_display),
            leave_out_level = leave_out_level,
            print_level = print_level
        )
    else
        settings = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=simulations),
            first_id_effects =first_id_effects,
            cov_effects = cov_effects,
            first_id_display = first_id_display,
            first_id_display_small = lowercase(first_id_display),
            second_id_display = second_id_display,
            second_id_display_small = lowercase(second_id_display),
            outcome_id_display= outcome_id_display,
            outcome_id_display_small = lowercase(outcome_id_display),
            leave_out_level = leave_out_level,
            print_level = print_level
        )
    end

    println("Running KSS correction with the following options")
    println("Leave Out Strategy : Leave $(settings.leave_out_level) out")
    if algorithm == "jla"
        algo = "JLA"
        simul = settings.leverage_algorithm.num_simulations == 0 ? 200 : settings.leverage_algorithm.num_simulations
        println("Algorithm for Computation of Statistical Leverages : $(algo) with $(simul) simulations")
    else
        algo = "Exact"
        println("Algorithm for Computation of Statistical Leverages : $(algo) \n")
    end

    first_id_old = first_id
    second_id_old = second_id
    y_untransformed = y

    @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y, first_id, second_id, settings, controls)

    #Residualize outcome variable
    if controls != nothing
        println("\nPartialling out controls...")
        NT = size(y,1)
        J = maximum(second_id)
        N = maximum(first_id)
        K = size(controls,2)
        nparameters = N + J + K

        D = sparse(collect(1:NT),first_id,1)
        F = sparse(collect(1:NT),second_id,1)
        S= sparse(1.0I, J-1, J-1)
        S=vcat(S,sparse(-zeros(1,J-1)))
        X = hcat(D, -F*S, controls)

        #My best shot is to wrap AMG as LinearOperator
        buff = zeros(size(X,2))
        xx = X'*X
        P = AmgOperator(ruge_stuben(xx),buff)
        xy=X'*y
        beta, stats = Krylov.cg(xx,[xy...]; M = P , rtol = 1e-6, itmax = 300)

        y=y-X[:,N+J:end]*beta[N+J:end]
        controls = nothing
        println("Partialling out completed.")
    end

    # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y,first_id,second_id,controls,settings)
    @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov, y, X, sigma_i = leave_out_estimation(y,first_id,second_id,controls,settings)

    Z_lincom = nothing
    if  parsed_args["do_lincom"]
        #Construct Transform
        F = sparse(collect(1:length(second_id)),second_id,1)
        J = size(F,2)
        S = sparse(1.0I, J-1, J-1)
        S = vcat(S,sparse(-zeros(1,J-1)))
        Transform = hcat(spzeros(length(second_id),maximum(first_id)), -F*S )

        #Construct Z_lincom
        if lincom_covariates == []
            println("\n User asked for lincom but no covariates were specified. This step will not be performed.")
        else
            if typeof(lincom_covariates[1]) == String
                #Build lincom controls matrix
                lincom_labels = []
                push!(lincom_labels,String.(lincom_covariates[1]))

                Z_lincom = data[obs,lincom_covariates[1]]
                if length(lincom_covariates)>=2
                    for k=2:length(lincom_covariates)
                        hcat(Z_lincom, data[obs,lincom_covariates[k]])
                        push!(lincom_labels,String.(lincom_covariates[k]))
                    end
                end
            else
                println("WARNING: Elements of lincom covariates are not defined correctly. No inference will be performed.")
            end
        end
    end

    if Z_lincom != nothing
        #Collapse and reweight to person-year observations
        match_id = compute_matchid(second_id, first_id)
        Z_lincom_col = ones(size(Z_lincom,1),1)
        for i = 1:size(Z_lincom,2)
            Z_lincom_col = hcat(Z_lincom_col,(transform(groupby(DataFrame(z = Z_lincom[:,i], match_id = match_id), :match_id), :z => mean  => :z_py).z_py))
        end

        @unpack test_statistic, linear_combination , SE_linear_combination_KSS, SE_naive = lincom_KSS(y,X, Z_lincom_col, Transform, sigma_i; lincom_labels)
    end

    if parsed_args["write_detailed_csv"]

        y_output = y
        first_id_output = first_id
        second_id_output = second_id

        max_length = length(obs)

        if leave_out_level == "match"
            #We need to take the output to the person-year space
            match_id = compute_matchid(second_id, first_id)
            weights = accumarray(match_id, 1)

            Pii = [Pii...]
            Pii = vcat(fill.(Pii, weights)...)
            Bii_first = Bii_first == nothing ? nothing : vcat(fill.(Bii_first, weights)...)
            Bii_second = Bii_second == nothing ? nothing : vcat(fill.(Bii_second, weights)...)
            Bii_cov = Bii_cov == nothing ? nothing : vcat(fill.(Bii_cov, weights)...)
        end

        #Make sure that we can put them in a DataFrame
        obs = [obs...]
        first_id_old = [first_id_old[obs]...]
        first_id_output = [first_id_output...]
        second_id_old = [second_id_old[obs]...]
        second_id_output = [second_id_output...]
        y_output = [y_untransformed[obs]...]
        Dalpha = [Dalpha...]
        Fpsi = [Fpsi...]
        Pii = [Pii...]

        #todo rename the DataFrame arguments
        output = DataFrame(observation = obs , first_id_old = first_id_old, first_id = first_id_output ,
                           second_id_old = second_id_old, second_id = second_id_output, y = y_output,
                           Dalpha = Dalpha , Fpsi = Fpsi, Pii = Pii,
                           Bii_first = Bii_first === nothing ? missings(max_length) : Bii_first,
                           Bii_second = Bii_second === nothing ? missings(max_length) :  Bii_second,
                           Bii_cov = Bii_cov === nothing ? missings(max_length) : Bii_cov
                        )
        output_path = parsed_args["detailed_csv_path"]
        #TODO not sure about the paths
        if isabspath(output_path) == false
            output_path = joinpath(pwd(), output_path)
        end
        CSV.write(output_path,output)
    end

    if parsed_args["write_results"]
        y_py = y_untransformed[obs]
        var_den = var(y_py)
        movers , T = compute_movers(first_id, second_id)
        num_movers = length(unique(movers .* first_id)) - 1

        output_template = """
            Number of observations (Leave Out Sample): $(length(obs)) \n
            Number of $(first_id_display_small)s: $(maximum(first_id)) \n
            Number of $(second_id_display_small)s: $(maximum(second_id)) \n
            Number of movers: $(num_movers) \n
            Mean Outcome: $(mean(y_py)) \n
            Variance Outcome: $(var(y_py)) \n
            Bias Corrected Variance of $(second_id_display) Effects: $θ_second \n
        """

        if parsed_args["write_results"]
            try
                output_path = parsed_args["results_path"]
                if isabspath(output_path) == false
                    output_path = joinpath(pwd(), output_path)
                end
                open(output_path, "w") do io
                   write(io, output_template)
                   if first_id_effects == 1
                        write(io, "    Bias Corrected Variance of  $(first_id_display) Effects:  $θ_first \n\n")
                   end

                   if cov_effects == 1
                        write(io, "    Bias Corrected Covariance of  $(first_id_display)-$(second_id_display) Effects: $θCOV \n\n")
                   end

                   if (first_id_effects == 1 ) && (cov_effects == 1 )
                        θ_first = abs(θ_first) #This is only required for the very small dataset that does the precompilation
                        θ_second = abs(θ_second) #This is only required for the very small dataset that does the precompilation
                        corr = θCOV/(sqrt(θ_first)*sqrt(θ_second))
                        r2 = (θ_second+2*θCOV+θ_first)/var_den
                        write(io, "    Bias Corrected Correlation of  $(first_id_display)-$(second_id_display) Effects: $corr \n\n")
                        write(io, "    Bias Corrected Fraction of Variance explained by  $(first_id_display)-$(second_id_display) Effects: $r2 \n\n")
                   end

                    if Z_lincom != nothing
                        #Write the output of inference
                        r = size(Z_lincom_col,2)
                        write(io,"    Results of High Dimensional Lincom \n\n")
                        if lincom_labels == nothing
                            for q=2:r
                                if q <= r
                                    ncol = q-1
                                    output_inference = """
                                        Coefficient of Column $(ncol): $(linear_combination[q]) \n
                                        Traditional HC Standard Error of Column $(ncol): $(SE_naive[q]) \n
                                        KSS Standard Error of Column $(ncol): $(SE_linear_combination_KSS[q]) \n
                                        T-Statistic of Column $(ncol): $(test_statistic[q]) \n
                                    """
                                    write(io,output_inference)
                                end
                            end
                        else
                            for q=2:r
                                tell_me = lincom_labels[q-1]
                                output_inference = """
                                    Coefficient on $(tell_me): $(linear_combination[q-1]) \n
                                    Traditional HC Standard Error on $(tell_me): $(SE_naive[q-1]) \n
                                    KSS Standard Error on $(tell_me): $(SE_linear_combination_KSS[q-1]) \n
                                    T-Statistic on $(tell_me): $(test_statistic[q-1]) \n
                                """
                                write(io,output_inference)
                            end
                        end
                    end
                 end
            catch e
                println("Error writing output to $(output_path)")
                display(e)
           end
        end
    end

end

end


text:
#!/usr/bin/env julia
#Start Test Script
using Test

#Test Case
@testset "Expand" begin
    tst = "a{name}b{number}c.ext"
    #test that we safely ignore unused key's/wildcards and get the right length
    #and also the correct result!
    semantics = Dict(   "name"    => [ "great", "okay", "bad" ],
                        "number"  => [ "1", "2", "3" ],
                        "slumber" => [ "z", "zz","zzz" ] )
    exper = [ i for i in Expand(tst,semantics) ]
    @test all( exper .== ["agreatb1c.ext", "aokayb1c.ext", "abadb1c.ext", "agreatb2c.ext", "aokayb2c.ext", "abadb2c.ext", "agreatb3c.ext", "aokayb3c.ext", "abadb3c.ext"] )
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
# # Log-linear Language Model

using Knet: Param, @diff, grad, params, value, KnetArray, logp, nll, progress
using Random: shuffle!

# We are using data from the Penn Treebank, which is already converted into an easy-to-use
# format with "<unk>" symbols. If we were using other data we would have to do
# pre-processing and consider how to choose unknown words, etc.

readdata(file)=[ w2i.(split(line)) for line in eachline(file) ]

# Before reading the training data, we initalize the word->id dictionary with sentence
# separator "<s>" and unknown word "<unk>" symbols.

wdict = Dict()
w2i(x) = get!(wdict, x, 1+length(wdict)) # insert key with value len+1 if not found
S = w2i("<s>")
UNK = w2i("<unk>")

# Load the training data and peek at the first instance:

trn = readdata("../data/ptb/train.txt")
first(trn)'

# Here is the reconstructed string of the first sentence:

wstring = Array{String}(undef,length(wdict))
for (str,id) in wdict; wstring[id] = str; end
i2w(i) = wstring[i]
join(i2w.(rand(trn)), " ")

# Before reading the dev/test data, we change the word->id function to return UNK for
# unknown words.

w2i(x) = get(wdict, x, UNK)     # return UNK if x is not found

# Load the dev/test data and print the number of instances:

dev = readdata("../data/ptb/valid.txt")
tst = readdata("../data/ptb/test.txt")
length.((trn, dev, tst))

# Use KnetArray to initialize parameters on GPU, use Array to initialize on CPU:

##param(dims...) = Param(Array(0.01f0 * randn(Float32, dims...)))
param(dims...) = Param(KnetArray(0.01f0 * randn(Float32, dims...)))

# Initialize the parameters of the loglin-lm model as global variables W and b.

nwords = length(wdict)
N = 2  # The length of the n-gram
W = [ param(nwords, nwords) for i in 1:N ]
b = param(nwords)

# Here is the loss function for a whole sentence:

function loss(sent)
    slen = length(sent)
    input = [ repeat([S],N); sent ]
    scores = b
    for i in 1:N
        scores = scores .+ W[i][:,input[i:i+slen]] # @size scores (V,slen+1)
    end
    nll(scores, [sent; [S]])
end

# Here is the SGD training loop:

function train(data=trn; nepochs = 10, lr = 0.01f0)
    for epoch in 1:nepochs
        shuffle!(data)
        lastloss = 0
        for s in progress(x->lastloss, data)
            ∇loss = @diff loss(s)
            lastloss = value(∇loss)
            for p in params(∇loss)
                p .= p - lr * grad(∇loss, p)
            end
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

using FEMBeam, Test
using FEMBeam: get_beam_stiffness_matrix_2d, get_beam_forces_vector_2d,
               get_beam_mass_matrix_2d

beam_element = Element(Seg2, [1, 2])
problem = Problem(Beam2D, "test problem", 3)
add_elements!(problem, [beam_element])

@test get_unknown_field_name(problem) == "displacement"
@test assemble!(problem, 0.0) == nothing

## Beam 1

X1 = [0.0, 0.0]
X2 = [0.0, 6.5]
E = 210.0e9
I1 = 50.8e-3*101.6e-3^3/12
A = 50.8e-3*101.6e-3
rho = 7800.0
qt = -500.0
qn = 0.0
fp = zeros(6)

k = get_beam_stiffness_matrix_2d(X1, X2, E, I1, A)
f = get_beam_forces_vector_2d(X1, X2, qt, qn, fp)
m = get_beam_mass_matrix_2d(X1, X2, A, rho)

m_expected = [
  97.1943   0.0     -89.0948   33.6442   0.0      52.6469
   0.0     87.2256    0.0       0.0     43.6128    0.0
 -89.0948   0.0     105.294   -52.6469   0.0     -78.9703
  33.6442   0.0     -52.6469   97.1943   0.0      89.0948
   0.0     43.6128    0.0       0.0     87.2256    0.0
  52.6469   0.0     -78.9703   89.0948   0.0     105.294]

k_expected = [
40740.3         0.0            -1.32406e5  -40740.3         0.0            -1.32406e5
    0.0         1.66749e8       0.0             0.0        -1.66749e8       0.0
   -1.32406e5   0.0        573759.0             1.32406e5   0.0             2.8688e5
-40740.3         0.0             1.32406e5   40740.3         0.0             1.32406e5
    0.0        -1.66749e8       0.0             0.0         1.66749e8       0.0
   -1.32406e5   0.0             2.8688e5        1.32406e5   0.0        573759.0]

f_expected = [1625.0, 0.0, -1760.42, 1625.0, 0.0, 1760.42]

@test isapprox(m, m_expected; rtol=1.0e-5)
@test isapprox(k, k_expected; rtol=1.0e-5)
@test isapprox(f, f_expected; rtol=1.0e-5)

## Beam 2

X1 = [0.0, 6.5]
X2 = [8.0, 6.5]
E = 250.0e9
I1 = 75.0e-3*100.0e-3^3/12
A = 75.0e-3*100.0e-3
rho = 10000.0
qt = -750.0
qn = 0.0
fp = zeros(6)

m = get_beam_mass_matrix_2d(X1, X2, A, rho)
k = get_beam_stiffness_matrix_2d(X1, X2, E, I1, A)
f = get_beam_forces_vector_2d(X1, X2, qt, qn, fp)

m_expected = [
200.0     0.0        0.0    100.0     0.0        0.0
  0.0   222.857    251.429    0.0    77.1429  -148.571
  0.0   251.429    365.714    0.0   148.571   -274.286
100.0     0.0        0.0    200.0     0.0        0.0
  0.0    77.1429   148.571    0.0   222.857   -251.429
  0.0  -148.571   -274.286    0.0  -251.429    365.714]

k_expected = [
  2.34375e8       0.0        0.0  -2.34375e8        0.0        0.0
  0.0         36621.1   146484.0   0.0         -36621.1   146484.0
  0.0        146484.0   781250.0   0.0        -146484.0   390625.0
 -2.34375e8       0.0        0.0   2.34375e8        0.0        0.0
  0.0        -36621.1  -146484.0   0.0          36621.1  -146484.0
  0.0        146484.0   390625.0   0.0        -146484.0   781250.0]

f_expected = [0.0, -3000.0, -4000.0, 0.0, -3000.0, 4000.0]

@test isapprox(m, m_expected; rtol=1.0e-5)
@test isapprox(k, k_expected; rtol=1.0e-5)
@test isapprox(f, f_expected; rtol=1.0e-5)

## Beam 3

X1 = [8.0, 6.5]
X2 = [8.0, 0.0]
E = 160.0e9
I1 = 50.0e-3*50.0e-3^3/12
A = 50.0e-3*50.0e-3
rho = 6000.0
qt = 0.0
qn = 0.0
fp = zeros(6)

m = get_beam_mass_matrix_2d(X1, X2, A, rho)
k = get_beam_stiffness_matrix_2d(X1, X2, E, I1, A)
f = get_beam_forces_vector_2d(X1, X2, qt, qn, fp)

m_expected = [
36.2143   0.0    33.1964   12.5357   0.0   -19.6161
 0.0     32.5     0.0       0.0     16.25    0.0
33.1964   0.0    39.2321   19.6161   0.0   -29.4241
12.5357   0.0    19.6161   36.2143   0.0   -33.1964
 0.0     16.25    0.0       0.0     32.5     0.0
-19.6161   0.0   -29.4241  -33.1964   0.0    39.2321]

k_expected = [
  3641.33   0.0         11834.3   -3641.33   0.0         11834.3
     0.0    6.15385e7       0.0       0.0   -6.15385e7       0.0
 11834.3    0.0         51282.1  -11834.3    0.0         25641.0
 -3641.33   0.0        -11834.3    3641.33   0.0        -11834.3
     0.0   -6.15385e7       0.0       0.0    6.15385e7       0.0
 11834.3    0.0         25641.0  -11834.3    0.0         51282.1]

f_expected = zeros(6)

@test isapprox(m, m_expected; rtol=1.0e-5)
@test isapprox(k, k_expected; rtol=1.0e-5)
@test isapprox(f, f_expected; rtol=1.0e-5)


text:
# 3D
function LinearAlgebra.eigvals(R::Rotation{3})
    θ = rotation_angle(R)
    return SVector(exp(-θ*im), exp(θ*im), 1)
end

function LinearAlgebra.eigvecs(R::Rotation{3})
    n3 = normalize(rotation_axis(R))
    n1 = normalize(perpendicular_vector(n3))
    n2 = normalize(n3×n1)
    v1 = normalize(n1 + n2*im)
    v2 = normalize(n1*im + n2)
    v3 = n3
    return hcat(v1,v2,v3)
end

function LinearAlgebra.eigen(R::Rotation{3})
    return eigen(AngleAxis(R))
end

function LinearAlgebra.eigen(R::AngleAxis)
    λs = eigvals(R)
    vs = eigvecs(R)
    return Eigen(λs, vs)
end

# 2D
function LinearAlgebra.eigvals(R::Rotation{2})
    θ = rotation_angle(R)
    return SVector(exp(-θ*im), exp(θ*im))
end

function LinearAlgebra.eigvecs(R::Rotation{2})
    return @SMatrix [1/√2 1/√2;im/√2 -im/√2]
end

function LinearAlgebra.eigen(R::Rotation{2})
    λs = eigvals(R)
    vs = eigvecs(R)
    return Eigen(λs, vs)
end


text:
using CUDAdrv

using Test

dev = CuDevice(0)
ctx = CuContext(dev)

md = CuModuleFile(joinpath(@__DIR__, "vadd.ptx"))
vadd = CuFunction(md, "kernel_vadd")

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = Mem.upload(a)
d_b = Mem.upload(b)
d_c = Mem.alloc(c)

len = prod(dims)
cudacall(vadd, Tuple{CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}}, d_a, d_b, d_c; threads=len)

Mem.download!(c, d_c)
@test a+b ≈ c

destroy!(ctx)


text:
JSONRPC.parse_params(::Type{Val{Symbol("textDocument/hover")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/hover")},TextDocumentPositionParams}, server)
    if !haskey(server.documents, URI2(r.params.textDocument.uri))
        send(JSONRPC.Response(r.id, CancelParams(r.id)), server)
        return
    end
    doc = server.documents[URI2(r.params.textDocument.uri)]
    x = get_expr(getcst(doc), get_offset(doc, r.params.position), 0, true)
    documentation = get_hover(x, "", server)

    send(JSONRPC.Response(r.id, Hover(MarkupContent(documentation), missing)), server)
end

function get_hover(x, documentation, server) documentation end

function get_hover(x::EXPR, documentation, server)
    if parentof(x) isa EXPR  && (kindof(x) === CSTParser.Tokens.END || kindof(x) === CSTParser.Tokens.RPAREN || kindof(x) === CSTParser.Tokens.RBRACE || kindof(x) === CSTParser.Tokens.RSQUARE)
        documentation = "Closes `$(typof(parentof(x)))` expression."
    elseif CSTParser.isidentifier(x) && StaticLint.hasref(x)
        if refof(x) isa StaticLint.Binding
            documentation = get_hover(refof(x), documentation, server)
        elseif refof(x) isa SymbolServer.SymStore
            documentation = string(documentation, refof(x).doc)
        end
    end
    return documentation
end

function get_hover(b::StaticLint.Binding, documentation, server)
    if b.val isa EXPR
        if CSTParser.defines_function(b.val)
            while true
                if b.val isa EXPR
                    if parentof(b.val) isa EXPR && typof(parentof(b.val)) === CSTParser.MacroCall && length(parentof(b.val).args) == 3 && typof(parentof(b.val).args[1]) === CSTParser.GlobalRefDoc && CSTParser.isstring(parentof(b.val).args[2])
                        # Binding has preceding docs so use them..
                        documentation = string(documentation, Expr(parentof(b.val).args[2]))
                    elseif CSTParser.defines_function(b.val)
                        documentation = string(documentation, "```julia\n", Expr(CSTParser.get_sig(b.val)), "\n```\n")
                    elseif CSTParser.defines_datatype(b.val)
                        documentation = string(documentation, "```julia\n", Expr(b.val), "\n```\n")
                    end
                elseif b.val isa SymbolServer.SymStore
                    documentation = string(documentation, b.val.doc)
                else
                    break
                end
                if b.prev isa StaticLint.Binding && b.prev != b && (b.prev.type == StaticLint.CoreTypes.Function || b.prev.type == StaticLint.CoreTypes.DataType || b.prev.val isa Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore})
                    b = b.prev
                else
                    break
                end
            end
        else
            documentation = string(documentation, "```julia\n", Expr(b.val), "\n```\n")
        end
    elseif b.val isa SymbolServer.SymStore
        documentation = string(documentation, b.val.doc)
    elseif b.val isa StaticLint.Binding
        documentation = get_hover(b.val, documentation, server)
    end
    return documentation
end


text:
# a simple LP example
#
#   maximize x + y
#
#   s.t. 50 x + 24 y <= 2400
#        30 x + 33 y <= 2100
#        x >= 45, y >= 5
#
#   solution: x = 45, y = 6.25, objv = 51.25

using Gurobi, Compat.Test, Compat.GC

@testset "LP 01b" begin

    env = Gurobi.Env()
    setparam!(env, "OutputFlag", 0)
    setparams!(env, Method=2)  # using barrier method

    method = getparam(env, "Method")

    model = Gurobi.Model(env, "lp_01", :maximize)

    # add variables
    add_cvars!(model, [1., 1.], [45., 5.], Inf)
    update_model!(model)

    # add constraints
    add_constrs!(model, Cint[1, 3], Cint[1, 2, 1, 2],
        [50., 24., 30., 33.], '<', [2400., 2100.])
    update_model!(model)

    # perform optimization
    optimize(model)

    # show results
    info = get_optiminfo(model)

    sol = get_solution(model)
    @test sol == [45, 6.25]

    objv = get_objval(model)
    @test objv == 51.25

    GC.gc()  # test finalizers
end


text:
function displace(x_prim::Float64, tau_x::Float64, lb_x::Float64, ub_x::Float64)

	delta_x::Float64 = tau_x * randn()
	x_bis::Float64 = x_prim + delta_x
	while !( lb_x < x_bis < ub_x )
		if x_bis < lb_x
			x_bis = x_bis + 2.0 * ( lb_x - x_bis )
		end
		if x_bis > ub_x
			x_bis = x_bis - 2.0 * ( x_bis - ub_x )
		end
	end

	return x_bis
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module BlapackPerf
import Perftests: @perf, meta

sizes = [("tiny", 1), ("small", 4), ("medium", 6), ("large", 8), ("huge", 10)]
for (size_str, exp) in sizes
    # The length of vectors, size of matrices, etc...
    n = 2^exp

    # BLAS level 1 tests
    a = rand(n)
    b = rand(n)
    @perf Base.dot(a,b) meta("dot", size_str, "Dot product with lengths 2^$exp")
    a = rand()
    x = rand(n)
    y = zeros(n)
    @perf Base.axpy!(a,x,y) meta("axpy", size_str, "Vector add-multiply with lengths 2^$exp")

    # BLAS level 2 tests
    A = rand(n,n)
    x = rand(n)
    @perf A*x meta("gemv", size_str, "Matrix-vector multiply with lengths 2^$exp")

    # BLAS level 3 tests
    B = similar(A)
    @perf A_mul_B!(B, A, A) meta("gemm", size_str, "Matrix-matrix multiply with lengths 2^$exp")

    # LAPACK eig tests
    AA = A + A'
    C = rand(n,n) + im*rand(n,n)
    @perf eig(A) meta("realeig", size_str, "Real matrix eigenfactorization with lengths 2^$exp")
    @perf eig(A + A') meta("symeig", size_str, "Symmetric matrix eigenfactorization with lengths 2^$exp")
    @perf eig(C + C') meta("hermeig", size_str, "Hermitian matrix eigenfactorization with lengths 2^$exp")

    # LAPACK factorization tests
    @perf svdfact(A) meta("svdfact", size_str, "Singular Value Decomposition with lengths 2^$exp")
    @perf schurfact(A) meta("schurfact", size_str, "Schur factorization with lengths 2^$exp")
    @perf cholfact(A'*A) meta("cholfact", size_str, "Cholesky factorization with lengths 2^$exp")
    @perf qrfact(A) meta("qrfact", size_str, "QR factorizaiton with lengths 2^$exp")
    @perf lufact(A) meta("lufact", size_str, "LU factorizaiton with lengths 2^$exp")
end

end # module


text:
__precompile__()

module RigidBodySim

include("util.jl")
include("core.jl")
include("control.jl")
include("visualization.jl")

@reexport using .Core
@reexport using .Control
@reexport using .Visualization

# Select DifferentialEquations exports
using DiffEqBase: ODEProblem, init, solve!, solve, CallbackSet
export ODEProblem, init, solve!, solve, CallbackSet

using DiffEqCallbacks: PeriodicCallback
export PeriodicCallback

using OrdinaryDiffEq: Tsit5, Vern7, RK4
export Tsit5, Vern7, RK4

end # module


text:
export calcFactorResidual

function approxConvBelief(dfg::AbstractDFG,
                          fc::DFGFactor,
                          target::Symbol,
                          measurement::AbstractVector=Tuple[];
                          solveKey::Symbol=:default,
                          N::Int=length(measurement),
                          skipSolve::Bool=false )
  #
  v1 = getVariable(dfg, target)
  N = N == 0 ? getNumPts(v1, solveKey=solveKey) : N
  # points and infoPerCoord
  pts, ipc = evalFactor(dfg, fc, v1.label, measurement, solveKey=solveKey, N=N, skipSolve=skipSolve)

  len = length(ipc)
  mask = 1e-14 .< abs.(ipc)
  partl = collect(1:len)[ mask ]

  # is the convolution infoPerCoord full or partial
  if sum(mask) == len
    # not partial
    return manikde!(getManifold(getVariable(dfg, target)), pts, partial=nothing)
  else
    # is partial
    return manikde!(getManifold(getVariable(dfg, target)), pts, partial=partl)
  end
end

approxConv(w...;kw...) = getPoints( approxConvBelief(w...;kw...), false)

"""
    $SIGNATURES

Calculate the sequential series of convolutions in order as listed by `fctLabels`, and starting from the
value already contained in the first variable.

Notes
- `target` must be a variable.
- The ultimate `target` variable must be given to allow path discovery through n-ary factors.
- Fresh starting point will be used if first element in `fctLabels` is a unary `<:AbstractPrior`.
- This function will not change any values in `dfg`, and might have slightly less speed performance to meet this requirement.
- pass in `tfg` to get a recoverable result of all convolutions in the chain.
- `setPPE` and `setPPEmethod` can be used to store PPE information in temporary `tfg`

DevNotes
- TODO strong requirement that this function is super efficient on single factor/variable case!
- FIXME must consolidate with `accumulateFactorMeans`
- TODO `solveKey` not fully wired up everywhere yet
  - tfg gets all the solveKeys inside the source `dfg` variables
- TODO add a approxConv on PPE option
  - Consolidate with [`accumulateFactorMeans`](@ref), `approxConvBinary`

Related

[`approxDeconv`](@ref), `LightDFG.findShortestPathDijkstra`, [`evalFactor`](@ref)
"""
function approxConvBelief(dfg::AbstractDFG,
                          from::Symbol,
                          target::Symbol,
                          measurement::AbstractVector=Tuple[];
                          solveKey::Symbol=:default,
                          N::Int = length(measurement),
                          tfg::AbstractDFG = initfg(),
                          setPPEmethod::Union{Nothing, Type{<:AbstractPointParametricEst}}=nothing,
                          setPPE::Bool= setPPEmethod !== nothing,
                          path::AbstractVector{Symbol}=Symbol[],
                          skipSolve::Bool=false  )
  #
  # @assert isVariable(dfg, target) "approxConv(dfg, from, target,...) where `target`=$target must be a variable in `dfg`"

  if from in ls(dfg, target)
    # direct request
    # TODO avoid this allocation for direct cases ( dfg, :x1x2f1, :x2[/:x1] )
    path = Symbol[from; target]
    varLbls = Symbol[target;]
  else
    # must first discover shortest factor path in dfg
    # TODO DFG only supports LightDFG.findShortestPathDijkstra at the time of writing (DFG v0.10.9)
    path = 0 == length(path) ? findShortestPathDijkstra(dfg, from, target) : path
    @assert path[1] == from "sanity check failing for shortest path function"

    # list of variables
    fctMsk = isFactor.(dfg, path)
    # which factors in the path
    fctLbls = path[fctMsk]
    # must still add
    varLbls =  union(lsf.(dfg, fctLbls)...)
    neMsk = exists.(tfg, varLbls) .|> x-> xor(x,true)
    # put the non-existing variables into the temporary graph `tfg`
    # bring all the solveKeys too
    addVariable!.(tfg, getVariable.(dfg, varLbls[neMsk]))
    # variables adjacent to the shortest path should be initialized from dfg
    setdiff(varLbls, path[xor.(fctMsk,true)]) .|> x->initManual!(tfg, x, getBelief(dfg, x))
  end

  # find/set the starting point
  idxS = 1
  pts = if varLbls[1] == from
    # starting from a variable
    getBelief(dfg, varLbls[1]) |> getPoints
  else
    # chain would start one later
    idxS += 1
    # get the factor
    fct0 = getFactor(dfg,from)
    # get the Matrix{<:Real} of projected points
    pts1Bel = approxConvBelief(dfg, fct0, path[2], measurement, solveKey=solveKey, N=N, skipSolve=skipSolve)
    if length(path) == 2
      return pts1Bel
    end
    getPoints(pts1Bel)
  end
  # didn't return early so shift focus to using `tfg` more intensely
  initManual!(tfg, varLbls[1], pts)
  # use in combination with setPPE and setPPEmethod keyword arguments
  ppemethod = setPPEmethod === nothing ? MeanMaxPPE : setPPEmethod
  !setPPE ? nothing : setPPE!(tfg, varLbls[1], solveKey, ppemethod)

  # do chain of convolutions
  for idx in idxS:length(path)
    if fctMsk[idx]
      # this is a factor path[idx]
      fct = getFactor(dfg, path[idx])
      addFactor!(tfg, fct)
      ptsBel = approxConvBelief(tfg, fct, path[idx+1], solveKey=solveKey, N=N, skipSolve=skipSolve)
      initManual!(tfg, path[idx+1], ptsBel)
      !setPPE ? nothing : setPPE!(tfg, path[idx+1], solveKey, ppemethod)
    end
  end

  # return target variable values
  return getBelief(tfg, target)
end

"""
    $(SIGNATURES)

Compute proposal belief on `vertid` through `fct` representing some constraint in factor graph.
Always full dimension variable node -- partial constraints will only influence subset of variable dimensions.
The remaining dimensions will keep pre-existing variable values.

Notes
- fulldim is true when "rank-deficient" -- TODO swap to false (or even float)
"""
function calcProposalBelief(dfg::AbstractDFG,
                            fct::DFGFactor,
                            target::Symbol,
                            measurement::AbstractVector=Tuple[];
                            N::Int=length(measurement),
                            solveKey::Symbol=:default,
                            dbg::Bool=false  )
  #
  # assuming it is properly initialized TODO
  proposal = approxConvBelief(dfg, fct, target, measurement, solveKey=solveKey, N=N)

  # return the proposal belief and inferdim, NOTE likely to be changed
  return proposal
end

# specifically the PartialPriorPassThrough dispatch
function calcProposalBelief(dfg::AbstractDFG,
                            fct::DFGFactor{<:CommonConvWrapper{<:PartialPriorPassThrough}},
                            target::Symbol,
                            measurement::AbstractVector=Tuple[];
                            N::Int=length(measurement),
                            solveKey::Symbol=:default,
                            dbg::Bool=false  )
  #

  # density passed through directly from PartialPriorPassThrough.Z
  fctFnc = getFactorType(fct)
  proposal = fctFnc.Z.heatmap.densityFnc

  # in case of partial, place the proposal into larger marginal/partial MKD
  proposal_ = if isPartial(fctFnc)
    # oldbel = getBelief(dfg, target, solveKey)
    varType = getVariableType(dfg, target)
    M = getManifold(varType)
    u0 = getPointIdentity(varType)
    # replace(oldbel, proposal)
    antimarginal(M,u0,proposal,Int[fctFnc.partial...])
  else
    proposal
  end

  # return the proposal belief and inferdim, NOTE likely to be changed
  return proposal_
end

"""
    $SIGNATURES

Compute the proposals of a destination vertex for each of `factors` and place the result
as belief estimates in both `dens` and `partials` respectively.

Notes
- TODO: also return if proposals were "dimension-deficient" (aka ~rank-deficient).
"""
function proposalbeliefs!(dfg::AbstractDFG,
                          destlbl::Symbol,
                          factors::AbstractVector{<:DFGFactor},
                          dens::AbstractVector{<:ManifoldKernelDensity},
                          measurement::AbstractVector=Tuple[];
                          solveKey::Symbol=:default,
                          N::Int=getSolverParams(dfg).N, #maximum([length(getPoints(getBelief(dfg, destlbl, solveKey))); getSolverParams(dfg).N]),
                          dbg::Bool=false  )
  #

  # populate the full and partial dim containers
  ipcs = Vector{Vector{Float64}}(undef, length(factors))

  vardim = getDimension(getVariable(dfg, destlbl))
  # get a proposal belief from each factor connected to destlbl
  for (count,fct) in enumerate(factors)
    ccwl = _getCCW(fct)
    # need way to convey partial information
    # determine if evaluation is "dimension-deficient" solvable dimension
    # FIXME, update to infoPerCoord
    fct_ipc = ones(vardim) # getFactorSolvableDim(dfg, fct, destlbl, solveKey)
    # convolve or passthrough to get a new proposal
    propBel_ = calcProposalBelief(dfg, fct, destlbl, measurement, N=N, dbg=dbg, solveKey=solveKey)
    # partial density
    propBel = if isPartial(ccwl)
      pardims = _getDimensionsPartial(ccwl)
      @assert [getFactorType(fct).partial...] == [pardims...] "partial dims error $(getFactorType(fct).partial) vs $pardims"
      AMP.marginal(propBel_, Int[pardims...])
    else
      propBel_
    end
    push!(dens, propBel)
    ipcs[count] = fct_ipc
  end
  # len = maximum(length.(ipcs))
  ipc = zeros(vardim)
  for _ipc in ipcs
    ipc .+= _ipc
  end

  ipc
end
# group partial dimension factors by selected dimensions -- i.e. [(1,)], [(1,2),(1,2)], [(2,);(2;)]

# WIP, see `_buildGraphByFactorAndTypes!` where pts are full MKD Beliefs, following #1351
# Legacy use in RoMEPlotting: plotFactor
# function approxConvBelief(fct::AbstractFactorRelative,
#                           varTypes::Union{<:Tuple,<:AbstractVector{<:InstanceType{T}}},
#                           mkds::Union{<:Tuple,<:AbstractVector{<:InstanceType{T}}};
#                           tfg::AbstractDFG=_buildGraphByFactorAndTypes!(fct,)
#                           ) where {T <: InferenceVariable}
#   #

# end

#


text:
function get_covariance(obj::Measurement)
    return jcall(obj, "getCovariance", RealMatrix, ())
end

function get_time(obj::Measurement)
    return jcall(obj, "getTime", jdouble, ())
end

function get_value(obj::Measurement)
    return jcall(obj, "getValue", RealVector, ())
end


text:
@add_problem constant function constant_Issue_166(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Issue #166
    α = Variable(5)
    fix!(α, ones(5, 1))

    # has const vexity, but not at the head
    c = (rand(5, 5) * α) * ones(1, 5)

    β = Variable(5)
    set_value!(β, ones(5))

    problem = minimize(sum(c * β), [β >= 0]; numeric_type = T)

    handle_problem!(problem)
    if test
        @test problem.optval ≈ evaluate(sum(c * β)) atol = atol rtol = rtol
        @test problem.optval ≈ 0.0 atol = atol rtol = rtol
        @test evaluate(β) ≈ zeros(5) atol = atol rtol = rtol
    end
end

@add_problem constant function constant_Issue_228(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2)
    y = Variable(2)
    fix!(x, [1 1]')
    prob = minimize(y' * (x + [2 2]'), [y >= 0]; numeric_type = T)

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.0 atol = atol rtol = rtol
    end

    prob = minimize(x' * y, [y >= 0]; numeric_type = T)

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.0 atol = atol rtol = rtol
    end
end

@add_problem constant function constant_Test_double_fix!(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    y = Variable()
    fix!(x, 1.0)
    prob = minimize(y * x, [y >= x, x >= 0.5]; numeric_type = T)

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 1.0 atol = atol rtol = rtol
    end

    fix!(x, 2.0)
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 4.0 atol = atol rtol = rtol
    end

    free!(x)
    fix!(y, 1.0)
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.5 atol = atol rtol = rtol
    end
end

@add_problem constant function constant_fix!_and_multiply(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    p = Variable()
    fix!(p, 1.0)
    x = Variable(2, 2)
    prob = minimize(tr(p * x), [x >= 1]; numeric_type = T)

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 2.0 atol = atol rtol = rtol
        @test evaluate(tr(p * x)) ≈ 2.0 atol = atol rtol = rtol
    end
end

@add_problem constant function constant_fix!_with_complex_numbers(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = ComplexVariable()
    fix!(x, 1.0 + im * 1.0)
    y = Variable()
    prob = minimize(
        real(x * y),
        [y >= 0.5, real(x) >= 0.5, imag(x) >= 0];
        numeric_type = T,
    )

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.5 atol = atol rtol = rtol
        @test evaluate(real(x * y)) ≈ 0.5 atol = atol rtol = rtol
        @test evaluate(y) ≈ 0.5 atol = atol rtol = rtol
    end

    free!(x)
    if test
        fix!(y)
    else # if we haven't solved the problem, we need to fix! it to a number
        fix!(y, 0.5)
    end
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.25 atol = atol rtol = rtol
        @test evaluate(real(x * y)) ≈ 0.25 atol = atol rtol = rtol
        @test real(evaluate(x)) ≈ 0.5 atol = atol rtol = rtol
        @test evaluate(y) ≈ 0.5 atol = atol rtol = rtol
    end

    if test
        @test_throws DimensionMismatch fix!(x, rand(2))
        @test_throws DimensionMismatch fix!(x, rand(2, 2))
    end
end

@add_problem constant function constant_fix!_with_vectors(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = ComplexVariable(5)
    fix!(x, ones(5) + im * ones(5))
    y = Variable()
    prob = minimize(
        real(y * sum(x)),
        [y >= 0.5, real(x) >= 0.5, imag(x) >= 0];
        numeric_type = T,
    )

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 2.5 atol = atol rtol = rtol
        @test evaluate(real(y * sum(x))) ≈ 2.5 atol = atol rtol = rtol
        @test evaluate(y) ≈ 0.5 atol = atol rtol = rtol
    end

    free!(x)
    if test
        fix!(y)
    else # if we haven't solved the problem, we need to fix! it to a number
        fix!(y, 0.5)
    end

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 1.25 atol = atol rtol = rtol
        @test evaluate(real(y * sum(x))) ≈ 1.25 atol = atol rtol = rtol
        @test real(evaluate(x)) ≈ 0.5 * ones(5) atol = atol rtol = rtol
        @test evaluate(y) ≈ 0.5 atol = atol rtol = rtol
    end

    if test
        @test_throws DimensionMismatch fix!(x, rand(5, 5))
        @test_throws DimensionMismatch fix!(x, rand(4))
    end
end


text:
using LinearAlgebra
using Espresso
import Espresso: @get, @get_or_create,
    genname, to_context, subs_bcast_with_dot,
    assign_chain

include("utils.jl")
include("rules.jl")
include("deriv.jl")
include("xdiff.jl")
include("xgrad0.jl")
include("helpers.jl")
include("debug.jl")


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf, JLD2, FileIO
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 4.1 Multidimensional Interpolation

p=50;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("1"); n=BigInt(300); dx=BigFloat("2")*L/(n-1);
x=convert(Array{BigFloat},(rand(n).-1/2)*1)
y=convert(Array{BigFloat},(rand(n).-1/2)*1)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end
V=TaylorBigF.calc_Vandermonde_2D(x,y,n,n)
@time Vinv=V\I

# @save string("ws",p,".jld2") p,L,n,dx,x,y,I,V,Vinv
# p=1000;setprecision(p); setrounding(BigFloat, RoundUp);
# @load string("ws",p,".jld2") p,L,n,dx,x,y,I,V,Vinv

# sos pigaine poli kalitera me tria matrixes gia x y xy me psiles dinameis, alla mono se eidikes sinartisis ph sinx + cosy
# V=calc_Vandermonde_2D_highO(x,y,n,BigInt(34))

function func(x,y) sin.(BigFloat("5")x).+cos.(exp.(BigFloat("2")y)) end

f=func(x,y)
# V=V[:,end:-1:1]
# a=V\f
a=Vinv*f

# plot(a)
# scatter3d(x,y,f); scatter3d!(x,y,V*a)
# maximum(abs.(f-V*a))
# xi=x.+BigFloat("0.01")
# yi=y.+BigFloat("0.01")
# Vi=calc_Vandermonde_2D_highO(xi,yi,n,BigInt(34))
xi=convert(Array{BigFloat},(rand(n).-1/2)*0.7)
yi=convert(Array{BigFloat},(rand(n).-1/2)*0.7)
Vi=TaylorBigF.calc_Vandermonde_2D(xi,yi,n,n)
fi_comp=Vi*a
fi_analyt=func(xi,yi)
diff_fi=maximum(abs.(fi_analyt-fi_comp));@sprintf("diff_fi=%1.3e",diff_fi)
scatter3d(xi,yi,fi_analyt,label="exact values"
    ,markershape=:circle,color=:black)
scatter3d!(xi,yi,fi_comp,label=string("approximated, p=",p)
    ,markershape=:xcross,legend=:topleft,color=:black)
# scatter3d!(xi,yi,fi_comp,label=string("approximated, p=",p)
#     ,markershape=:xcross,legend=:topleft,color=:black,xaxis="x")
# zlims!(-1.4,2.0)
# maximum(fi_comp)
# minimum(fi_comp)
cd(dirname(@__FILE__))
savefig(string("interp3Dp",string(p),".pdf"))
# mean(abs.(V))-mean(abs.(Vi))
# mean(abs.(x))-mean(abs.(xi))


text:
include("ndid_core_supp.jl")

function eco_CR!(du, u, p, t)
    @unpack In1, In2, rn, a1, a11, a23, a2, a42, a43, h42, h43, b1, b11, b2, f, e1, e2, e3, dD, df, dc, d1, d11, d2 = p

    du[1] = In1 - rn * u[1] -a1 * u[1] * u[2] / (b1 + u[1]) - a11 * u[1] * u[3] / (b11 + u[1]) + dD * u[5]
    du[2] = f * a1 * u[1] * u[2] / (b1 + u[1]) - a42 * u[2] * u[4] / (1. + a42*h42*u[2] + a43*h43*u[3]) - (d1 + e1) * u[2]
    du[3] = f * a11 * u[1] * u[3] / (b11 + u[1]) -  a43 * u[3] * u[4] / (1. + a42*h42*u[2] + a43*h43*u[3]) - (d11 + e1) * u[3]
    du[4] = f * a42 * u[2] * u[4] / (1. + a42*h42*u[2] + a43*h43*u[3]) +  f * a43 * u[3] * u[4] / (1. + a42*h42*u[2] + a43*h43*u[3]) - (d2 + e2) * u[4] - df * u[4]
    du[5] = d1 * u[2] + d11 * u[3] + d2 * u[4] -  (dD + e3) * u[5]

    du[6] = In1 - rn * u[6] -a1 * u[6] * u[7] / (b1 + u[6]) - a11 * u[6] * u[8] / (b11 + u[6]) + dD * u[10]
    du[7] = f * a1 * u[6] * u[7] / (b1 + u[6]) - a42 * u[7] * u[9] / (1. + a42*h42*u[7] + a43*h43*u[8]) - (d1 + e1) * u[7]
    du[8] = f * a11 * u[6] * u[8] / (b11 + u[6]) - a43 * u[8] * u[9] / (1. + a42*h42*u[7] + a43*h43*u[8]) - (d11 + e1) * u[8]
    du[9] = f * a42 * u[7] * u[9] / (1. + a42*h42*u[7] + a43*h43*u[8]) + f * a43 * u[8] * u[9] / (1. + a42*h42*u[7] + a43*h43*u[8]) - (d2 + e2) * u[9] - df * u[9]
    du[10] = d1 * u[7] + d11 * u[8] + d2 * u[9] -  (dD + e3) * u[10]

    du[11] = In1 - rn * u[11] -a1 * u[11] * u[12] / (b1 + u[11]) - a11 * u[11] * u[13] / (b11 + u[11]) + dD * u[15]
    du[12] = f * a1 * u[11] * u[12] / (b1 + u[11]) - a42 * u[12] * u[14] / (1. + a42*h42*u[12] + a43*h43*u[13]) - (d1 + e1) * u[12]
    du[13] = f * a11 * u[11] * u[13] / (b11 + u[11])  - a43 * u[13] * u[14] / (1. + a42*h42*u[12] + a43*h43*u[13]) - (d11 + e1) * u[13]
    du[14] = f * a42 * u[12] * u[14] / (1. + a42*h42*u[12] + a43*h43*u[13]) + f * a43 * u[13] * u[14] / (1. + a42*h42*u[12] + a43*h43*u[13]) - (d2 + e2) * u[14] - df * u[14]
    du[15] = d1 * u[12] + d11 * u[13] + d2 * u[14] -  (dD + e3) * u[15]

# division by 3 to acount for 3fold increase in water volume
    du[16] = In1 - rn * u[16] -a1 * u[16] * u[17] / (b1 + u[16]) - a11 * u[16] * u[18] / (b11 + u[16]) + dD * u[20]
    du[17] = f * a1 * u[16] * u[17] / (b1 + u[16]) - a42 * u[17] * u[19] / (1. + a42*h42*u[17] + a43*h43*u[18]) - (d1 + e1) * u[17]
    du[18] = f * a11 * u[16] * u[18] / (b11 + u[16]) - a43 * u[18] * u[19] / (1. + a42*h42*u[17] + a43*h43*u[18]) - (d11 + e1) * u[18]
    du[19] = f * a42 * u[17] * u[19] / (1. + a42*h42*u[17] + a43*h43*u[18]) + f * a43 * u[18] * u[19] / (1. + a42*h42*u[17] + a43*h43*u[18]) - (d2 + e2) * u[19] + df * (u[4] + u[9] + u[14])/3. - df * u[19]
    du[20] = d1 * u[17] + d11 * u[18] + d2 * u[19] -  (dD + e3) * u[20]

    du[21] = In1 - rn * u[21] -a1 * u[21] * u[22] / (b1 + u[21]) - a11 * u[21] * u[23] / (b11 + u[21]) + dD * u[25]
    du[22] = f * a1 * u[21] * u[22] / (b1 + u[21]) - a42 * u[22] * u[24] / (1. + a42*h42*u[22] + a43*h43*u[23]) - (d1 + e1) * u[22]
    du[23] = f * a11 * u[21] * u[23] / (b11 + u[21]) - a43 * u[23] * u[24] / (1. + a42*h42*u[22] + a43*h43*u[23]) - (d11 + e1) * u[23]
    du[24] = f * a42 * u[22] * u[24] / (1. + a42*h42*u[22] + a43*h43*u[23]) + f * a43 * u[23] * u[24] / (1. + a42*h42*u[22] + a43*h43*u[23]) - (d2 + e2) * u[24] + df*u[19]
    du[25] = d1 * u[22] + d11 * u[23] + d2 * u[24] -  (dD + e3) * u[25]

    return du
end

@with_kw mutable struct EcoPar
    In1 = 0.15
    In2 = 0.0005
    rn = 0.05
    f = 0.8
    e1 = 0.05
    e2 = 0.05
    e3 = 0.05
    a1 = 0.6
    a11 = 0.55
    a23 = 0.07
    a2 = 0.25
    h42=.0005
    h43=.005
    b1 = 0.040
    b11 = 0.043
    b2 = 0.040
    a42 = 0.350
    a43 = 0.05
    dD = 0.1
    dfd = 0.0
    dcd = 0.0
    d1 = 0.10
    d11 = 0.11
    df2 = 0.001
    df = 0.2
    dc = 0.01
    d2 = 0.001
end

# Setup a default set
param = EcoPar()

# Example of simulation
u0 = fill(0.3, 25)
tspan = (0.0, 10000.0)
prob = ODEProblem(eco_CR!, u0, tspan, param)
sol = solve(prob, KenCarp5(), reltol = 1e-16)
plot(sol.t, sol.u)

#C's
@show "nutrients"
@show sol[1, end]
@show sol[6, end]
@show sol[11, end]
@show sol[16, end]
@show sol[21, end]

#Rs
@show "resources"
@show sol[2, end]
@show sol[7, end]
@show sol[12, end]
@show sol[17, end]
@show sol[22, end]

#RInedibles
@show "resources inedible"
@show sol[3, end]
@show sol[8, end]
@show sol[13, end]
@show sol[18, end]
@show sol[23, end]
#Cs
@show "consumers"
@show sol[4, end]
@show sol[9, end]
@show sol[14, end]
@show sol[19, end]
@show sol[24, end]

#Ds
@show "detritus"
@show sol[5, end]
@show sol[10, end]
@show sol[15, end]
@show sol[20, end]
@show sol[25, end]

plot(sol.t, sol[23, :])

function bif_analysis(vals)
    tspan = (0.0, 10000.0)
    u0 = fill(0.03, 25)
    # inedibles non-zero

    pbif = EcoPar()

    apoints3 = []
    apoints23 = []
    apoints4 = []
    apoints5 = []
    apoints24 = []
    apoints25 = []
    apoints2 = []
    apoints17 = []
    apoints22 = []
    acv3 = []
    acv23 = []

    for df in vals
        pbif.df = df
        prob = ODEProblem(eco_CR!, u0, tspan, pbif)
        sol = solve(prob, KenCarp5(), reltol = 1e-16)
        #@show sol.t
        maxes3 = local_max(sol(5000:1:10000)[3, :])
        mins3 = local_min(sol(5000:1:10000)[3, :])

        maxes23 = local_max(sol(5000:1:10000)[23, :])
        mins23 = local_min(sol(5000:1:10000)[23, :])

        maxes2 = local_max(sol(5000:1:10000)[2, :])
        mins2 = local_min(sol(5000:1:10000)[2, :])

        maxes17 = local_max(sol(5000:1:10000)[17, :])
        mins17 = local_min(sol(5000:1:10000)[17, :])

        maxes22 = local_max(sol(5000:1:10000)[22, :])
        mins22 = local_min(sol(5000:1:10000)[22, :])

        maxes4 = local_max(sol(5000:1:10000)[4, :])
        mins4 = local_min(sol(5000:1:10000)[4, :])

        maxes5 = local_max(sol(5000:1:10000)[5, :])
        mins5 = local_min(sol(5000:1:10000)[5, :])

        maxes24 = local_max(sol(5000:1:10000)[24, :])
        mins24 = local_min(sol(5000:1:10000)[24, :])

        maxes25 = local_max(sol(5000:1:10000)[25, :])
        mins25 = local_min(sol(5000:1:10000)[25, :])

        SD23 = std(sol(5000:1:10000)[23, :])
        MN23 = mean(sol(5000:1:10000)[23, :])
        CV23 = SD23 / MN23

        SD3 = std(sol(5000:1:10000)[3, :])
        MN3 = mean(sol(5000:1:10000)[3, :])
        CV3 = SD3 / MN3

        for m in maxes3
            push!(apoints3, [df, m])
        end
        for l in mins3
            push!(apoints3, [df, l])
        end

        for m in maxes23
            push!(apoints23, [df, m])
        end
        for l in mins23
            push!(apoints23, [df, l])
        end

        for m in maxes22
            push!(apoints22, [df, m])
        end
        for l in mins22
            push!(apoints22, [df, l])
        end

        for m in maxes17
            push!(apoints17, [df, m])
        end
        for l in mins17
            push!(apoints17, [df, l])
        end

        for m in maxes2
            push!(apoints2, [df, m])
        end
        for l in mins2
            push!(apoints2, [df, l])
        end

        for m in maxes4
            push!(apoints4, [df, m])
        end
        for l in mins4
            push!(apoints4, [df, l])
        end

        for m in maxes24
            push!(apoints24, [df, m])
        end
        for l in mins24
            push!(apoints24, [df, l])
        end

        for m in maxes25
            push!(apoints25, [df, m])
        end
        for l in mins25
            push!(apoints25, [df, l])
        end

        for m in maxes5
        push!(apoints5, [df, m])
        end
        for l in mins5
        push!(apoints5, [df, l])
        end

        for m in CV23
            push!(acv23, [df, m])
        end
        for m in CV3
            push!(acv3, [df, m])
        end

        @show df
    end

    return VectorOfArray(apoints2),VectorOfArray(apoints17),VectorOfArray(apoints22)
end

# Fig.S.5.B.ii (model 2)

bif = bif_analysis(0.0:0.001:.10)

let
    figure()

    subplot(131)
    plot(bif[1][1, :], bif[1][2, :], "ko", markersize = 1.4,color = "black")
    xlabel("df", fontname = "Times New Roman", fontsize = 12)
    ylabel("Max & Min R")
    xlim(0.0,.05)
    ylim(0.0,0.2)

    subplot(132)
    plot(bif[2][1, :], bif[2][2, :], "ko", markersize = 1.4)
    xlabel("df", fontname = "Times New Roman", fontsize = 12)
    ylabel("Max & Min R")
    xlim(0.0,.05)
    ylim(0.0,0.2)

    subplot(133)
    plot(bif[3][1, :], bif[3][2, :], "ko", markersize = 1.4)
    xlabel("df", fontname = "Times New Roman", fontsize = 12)
    ylabel("Max & Min R")
    xlim(0.0,.05)
    ylim(0.0,0.2)

    tight_layout()

end

max_stars_repo_name: McCannLab/NDID

max_issues_repo_name: McCannLab/NDID

max_forks_repo_name: McCannLab/NDID


text:
using Base.BinaryPlatforms

"""
    check_binfmt_misc_loaded()

Check that the `binfmt_misc` kernel module is loaded and enabled.
"""
function check_binfmt_misc_loaded()
    # If we're not running on Linux, then clearly `binfmt_misc` is not available
    if !Sys.islinux()
        return false
    end

    # If the `binfmt_misc` directory doesn't exist, the kernel module is likely not installed
    if !isdir("/proc/sys/fs/binfmt_misc")
        return false
    end

    # If the `status` file does not exist, the kernel module may not be loaded, or the
    # special `binfmt_misc` filesystem may not be mounted.
    if !isfile("/proc/sys/fs/binfmt_misc/status")
        return false
    end

    # Finally, check that the module itself has not been disabled.
    return strip(String(read("/proc/sys/fs/binfmt_misc/status"))) == "enabled"
end

"""
    BinFmtRegistration

Provides a structured view of a `binfmt_misc` interpreter registration.  Note that only "magic"
matching rules are allowed, we do not support "extension" matching rules.
"""
struct BinFmtRegistration
    name::String
    interpreter::String
    flags::Vector{Symbol}
    offset::Int64
    magic::Vector{UInt8}
    mask::Vector{UInt8}

    function BinFmtRegistration(name::AbstractString,
                                interpreter::AbstractString,
                                flags::Union{AbstractString,Vector{Symbol}},
                                offset::Integer,
                                magic::Vector{UInt8},
                                mask::Union{Nothing,Vector{UInt8}} = nothing)
        # Default mask is all `0xff`.
        if mask === nothing
            mask = UInt8[0xff for _ in 1:length(magic)]
        end
        if isa(flags, AbstractString)
            flags = Symbol.(collect(flags))
        end
        return new(String(name), String(interpreter), sort(flags), Int64(offset), magic, mask)
    end
end

"""
    register_string(reg::BinFmtRegistration)

Constructs the string used to register a `binfmt_misc` registration with the `register`
file endpoint within `/proc/sys/fs/binfmt_misc/register`.  To actually register the
interpreter, use `write_binfmt_misc_registration()`.
"""
function register_string(reg::BinFmtRegistration)
    return string(
        ":",
        reg.name,
        ":",
        # We only support `magic` style registrations
        "M",
        ":",
        string(reg.offset),
        ":",
        # We need to actually emit double-escaped hex codes, since that's what `/register` expects.
        join([string("\\x", string(x, base=16, pad=2)) for x in reg.magic], ""),
        ":",
        join([string("\\x", string(x, base=16, pad=2)) for x in reg.mask], ""),
        ":",
        reg.interpreter,
        ":",
        join(String.(reg.flags), ""),
    )
end

macro check_specified(name)
    return quote
        if $(esc(name)) === nothing
            throw(ArgumentError($("Error, $(name) must be specified")))
        end
    end
end

"""
    BinFmtRegistration(file::String)

Reads a `binfmt_misc` registration in from disk, if it cannot be parsed (because it is
malformed, or uses unsupported features) it an `ArgumentError` will be thrown.
"""
function BinFmtRegistration(file::String)
    enabled = false
    interpreter = nothing
    flags = nothing
    offset = nothing
    magic = nothing
    mask = nothing
    for l in strip.(filter(!isempty, split(String(read(file)), "\n")))
        # Handle enabled/disabled line
        if l in ("enabled", "disabled")
            enabled = l == "enabled"
        elseif startswith(l, "interpreter ")
            interpreter = l[13:end]
        elseif startswith(l, "flags:")
            flags = l[8:end]
        elseif startswith(l, "offset ")
            offset = parse(Int64, l[8:end])
        elseif startswith(l, "magic ")
            magic = hex2bytes(l[7:end])
        elseif startswith(l, "mask ")
            mask = hex2bytes(l[6:end])
        else
            @warn("Unknown `binfmt_misc` configuration directive", line=l)
        end
    end

    # Ensure we are only dealing with properly fully-specified binfmt_misc registrations
    @check_specified interpreter
    @check_specified flags
    @check_specified offset
    @check_specified magic

    # If we found a disabled binfmt_misc registration, just ignore it
    if !enabled
        return nothing
    end

    return BinFmtRegistration(basename(file), interpreter, flags, offset, magic, mask)
end

function formats_match(a::BinFmtRegistration, b::BinFmtRegistration)
    return (a.magic .& a.mask) == (b.magic .& b.mask)
end

"""
    read_binfmt_misc_registrations()

Return a list of `BinFmtRegistration` objects, one per readable registration, as found
sitting in `/proc/sys/fs/binfmt_misc/*`.  Registrations that cannot be parsed are
silently ignored.
"""
function read_binfmt_misc_registrations()
    if !check_binfmt_misc_loaded()
        return String[]
    end

    registrations = BinFmtRegistration[]
    for f in readdir("/proc/sys/fs/binfmt_misc"; join=true)
        # Skip "special" files
        if basename(f) ∈ ("register", "status")
            continue
        end

        try
            reg = BinFmtRegistration(f)
            if reg !== nothing
                push!(registrations, reg)
            end
        catch e
            if isa(e, ArgumentError)
                continue
            end
            rethrow(e)
        end
    end
    return registrations
end

sudo_tee(f::Function, path::String) = open(f, Cmd([sudo_cmd()..., "tee", "-a", path]), write=true)

"""
    write_binfmt_misc_registration(reg::BinFmtRegistration)

Write a `binfmt_misc` registration out to the kernel's `register` file endpoint.
Requires `sudo` privileges.
"""
function write_binfmt_misc_registration!(reg::BinFmtRegistration)
    try
        sudo_tee("/proc/sys/fs/binfmt_misc/register") do io
            write(io, register_string(reg))
        end
    catch e
        @error("Unable to register binfmt_misc format", register_string=register_string(reg))
        rethrow(e)
    end
end

function clear_binfmt_misc_registrations!()
    sudo_tee("/proc/sys/fs/binfmt_misc/status") do io
        write(io, "-1")
    end
    return nothing
end

"""
    register_requested_formats(formats::Vector{BinFmtRegistration})

Given the list of `binfmt_misc` formats, check the currently-registered formats through
`read_binfmt_misc_registrations()`, check to see if any in `formats` are not yet
registered, and if they are not, call `write_binfmt_misc_registration!()` to register
it with an artifact-sourced `qemu-*-static` binary.
"""
function register_requested_formats!(formats::Vector{BinFmtRegistration}; verbose::Bool = false)
    # Do nothing if we're not asking for any formats.
    if isempty(formats)
        return nothing
    end

    # Read in the current binfmt_misc registrations:
    if !check_binfmt_misc_loaded()
        error("Cannot provide multiarch support if `binfmt_misc` not loaded!")
    end
    regs = read_binfmt_misc_registrations()

    # For each format, If there are no pre-existing registrations add it to `formats_to_register`
    formats_to_register = BinFmtRegistration[]
    for reg in formats
        if !any(formats_match.(Ref(reg), regs))
            push!(formats_to_register, BinFmtRegistration(
                reg.name,
                # We need to fetch our `multiarch-support` artifact, which has the necessary `qemu` executable.
                @artifact_str("multiarch-support/$(reg.name)-static"),
                reg.flags,
                reg.offset,
                reg.magic,
                reg.mask,
            ))
        end
    end

    # Notify the user if we have any formats to register, then register them.
    if !isempty(formats_to_register)
        format_names = sort([f.name for f in formats_to_register])
        msg = "Registering $(length(formats_to_register)) binfmt_misc entries, this may ask for your `sudo` password."
        if verbose
            @info(msg, formats=format_names)
        elseif (Sys.which("sudo") !== nothing) && (success(`sudo -k -n true`))
            # in this case, we know that `sudo` will not prompt the user for a password
        else
            @info(msg)
        end
        write_binfmt_misc_registration!.(formats_to_register)
    end
    return nothing
end

## binfmt_misc registration templates for various architectures.
## Note that these are true no matter the host architecture; e.g. these
## can just as easily point at `x86_64-qemu-aarch64-static` as `ppc64le-qemu-aarch64-static`.
## In fact, the interpreter path typically gets overwritten in `build_executor_command` anyway.
const qemu_x86_64 = BinFmtRegistration(
    "qemu-x86_64",
    "/usr/bin/qemu-x86_64-static",
    "OFC",
    0,
    UInt8[0x7f, 0x45, 0x4c, 0x46, 0x02, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x3e, 0x00],
    UInt8[0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe, 0xff, 0xff, 0xff],
)
const qemu_i386 = BinFmtRegistration(
    "qemu-i386",
    "/usr/bin/qemu-i386-static",
    "OFC",
    0,
    UInt8[0x7f, 0x45, 0x4c, 0x46, 0x02, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x03, 0x00],
    UInt8[0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe, 0xff, 0xff, 0xff],
)
const qemu_aarch64 = BinFmtRegistration(
    "qemu-aarch64",
    "/usr/bin/qemu-aarch64-static",
    "OFC",
    0,
    UInt8[0x7f, 0x45, 0x4c, 0x46, 0x02, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0xb7, 0x00],
    UInt8[0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe, 0xff, 0xff, 0xff],
)
const qemu_arm = BinFmtRegistration(
    "qemu-arm",
    "/usr/bin/qemu-arm-static",
    "OFC",
    0,
    UInt8[0x7f, 0x45, 0x4c, 0x46, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x28, 0x00],
    UInt8[0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe, 0xff, 0xff, 0xff],
)
const qemu_ppc64le = BinFmtRegistration(
    "qemu-ppc64le",
    "/usr/bin/qemu-ppc64le-static",
    "OFC",
    0,
    UInt8[0x7f, 0x45, 0x4c, 0x46, 0x02, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x15, 0x00],
    UInt8[0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfc, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe, 0xff, 0xff, 0x00],
)

const platform_qemu_registrations = Dict(
    # We register these `qemu-*-static` executables as capable of interpreting both `glibc` and `musl` platforms:
    Platform("x86_64", "linux"; libc="glibc") => qemu_x86_64,
    Platform("x86_64", "linux"; libc="musl") => qemu_x86_64,
    Platform("i686", "linux"; libc="glibc") => qemu_i386,
    Platform("i686", "linux"; libc="musl") => qemu_i386,
    Platform("aarch64", "linux"; libc="glibc") => qemu_aarch64,
    Platform("aarch64", "linux"; libc="musl") => qemu_aarch64,
    Platform("armv7l", "linux"; libc="glibc") => qemu_arm,
    Platform("armv7l", "linux"; libc="musl") => qemu_arm,
    Platform("ppc64le", "linux"; libc="glibc") => qemu_ppc64le,
    Platform("ppc64le", "linux"; libc="musl") => qemu_ppc64le,
)

# Define what is a natively-runnable
const host_arch = arch(HostPlatform())
function natively_runnable(p::Platform)
    if host_arch == "x86_64"
        return arch(p) ∈ ("x86_64", "i686")
    end
    if host_arch == "aarch64"
        return arch(p) ∈ ("aarch64", "armv7l")
    end
    return host_arch == arch(p)
end


text:
function benchmark_momentum(timestep::Vector;
    Nsim=1,
    tsim=1.0,
    tctrl=1.0,
    seed=0,
    ϵ=1.0e-14,
    control_amplitude=0.05)

    Nt = length(timestep)
    speed = zeros(Nt)
    plin = zeros(Nt)
    pang = zeros(Nt)
    for i = 1:Nt
        @show timestep[i]
        mech, storage, tcompute = astronaut_simulation(
            Nsim=Nsim,
            timestep=timestep[i],
            tsim=tsim,
            tctrl=tctrl,
            seed=seed,
            ϵ=ϵ,
            control_amplitude=control_amplitude)
        speed[i], plin[i], pang[i] = process_momentum_runs(mech, storage, tcompute)
    end
    return speed, plin, pang
end

function benchmark_energy(timestep::Vector;
    Nsim=1,
    tsim=2.0,
    tctrl=1.0,
    seed=0,
    ϵ=1.0e-14,
    control_amplitude=0.05)

    Nt = length(timestep)
    speed = zeros(Nt)
    energy = zeros(Nt)
    for i = 1:Nt
        @show timestep[i]
        mech, storage, tcompute = astronaut_simulation(
            Nsim=Nsim,
            timestep=timestep[i],
            tsim=tsim,
            tctrl=tctrl,
            seed=seed,
            ϵ=ϵ,
            control_amplitude=control_amplitude)
        speed[i], energy[i] = process_energy_runs(mech, storage, tcompute, tctrl)
    end
    return speed, energy
end


text:
"""
Remap a value from `(low1, high1)` to `(low2, high2)`.
"""
function remap(value, low1, high1, low2, high2)
    low2 + (value - low1) * (high2 - low2) / (high1 - low1)
end

remap(value, from, to) = remap(value, from..., to...)
remap(values::AbstractArray, to; from = extrema(values)) = remap.(values, Ref(from), Ref(to))

remap(low1, high1, low2, high2) = x -> remap(x, low1, high1, low2, high2)

not_implemented_for(T) = error("Not implemented for ", T)

"""
    @forward MyType.prop method1, method2, ...

Extend the provided methods by forwarding the property `prop` of `MyType` instances.
This will give, for a given `method`:
```julia
method(x::MyType, args...; kwargs...) = method(x.prop, args...; kwargs...)
```

"""
macro forward(ex, fs)
    T, prop = @match ex begin
        :($T.$prop) => (T, prop)
        _ => error("Invalid expression $ex, expected <Type>.<prop>")
    end

    fs = @match fs begin
        :(($(fs...),)) => fs
        :($mod.$method) => [fs]
        ::Symbol => [fs]
        _ => error("Expected a method or a tuple of methods, got $fs")
    end

    defs = map(fs) do f
        esc(:($f(x::$T, args...; kwargs...) = $f(x.$prop, args...; kwargs...)))
    end

    Expr(:block, defs...)
end


text:
mutable struct TBLogger <: AbstractLogger
    logdir::String
    file::Union{IOStream, Nothing}
    all_files::Dict{String, IOStream}
    global_step::Int
    step_increment::Int
    min_level::LogLevel
end

"""
    InitPolicy

An enum that describes how to behave if the folder that tensorboard should log
already exists.
"""
@enum InitPolicy tb_append=1 tb_overwrite=2 tb_increment=3 tb_read=4
export tb_append, tb_overwrite, tb_increment, tb_read

"""
    TBLogger(logdir[, tb_increment];
            time=time(),
            purge_step=nothing,
            step_increment=1,
            min_level=Logging.Info)

Creates a TensorBoardLogger in the folder `logdir`. The second (optional)
argument specifies the behaviour if the `logdir` already exhists: the default
choice `tb_increment` appends an increasing number 1,2... to `logdir`. Other
choices are `tb_overwrite`, which overwrites the previous folder, and `tb_append`.

If a `purge_step::Int` is passed, every step before `purge_step` will be ignored
by tensorboard (usefull in the case of restarting a crashed computation).

`min_level` specifies the minimum level of messages logged to
tensorboard.
"""
function TBLogger(logdir="tensorboard_logs/run", overwrite=tb_increment;
                  time=time(),
                  purge_step::Union{Int,Nothing}=nothing,
                  step_increment = 1,
                  min_level::LogLevel=Info)

    logdir = init_logdir(logdir, overwrite)

    if overwrite != tb_read
        fpath, evfile = create_eventfile(logdir, purge_step, time)
        all_files  = Dict(fpath => evfile)
    else
        evfile = nothing
        all_files = Dict{String, IOStream}()
    end

    start_step = something(purge_step, 0)

    TBLogger(logdir, evfile, all_files, start_step, step_increment, min_level)
end

"""
	init_logdir(logdir, [overwrite=tb_increment])

Creates a folder at path `logdir`. If the folder already exhists the behaviour
is determined by `overwrite`.

 - `overwrite=tb_increment` appends an increasing number 1,2... to logdir.
 - `overwrite=tb_overwrite` overwrites the folder, deleting it's content.
 - `overwrite=tb_append` appends to it's previous content.
 - `overwrite=tb_read` reads for deserialization only.
"""
function init_logdir(logdir, overwrite=tb_increment)
    if overwrite == tb_overwrite
        rm(logdir; force=true, recursive=true)
    elseif overwrite == tb_append
        # do nothing
    elseif overwrite == tb_read
        # logdir must exist
        ispath(logdir) || throw(ArgumentError("Logdir $logdir does not exist!"))
    elseif overwrite == tb_increment
        # find the next viable name
        if ispath(logdir)
            i=1
            while i<1000 # avoids an infinite loop
                logdir_new = logdir*"_$i"
                #println("trying $logdir_new")
                !ispath(logdir_new) && break
                i+=1
            end
            logdir = logdir*"_$i"
        end
    end
    mkpath(logdir)

    return realpath(logdir)
end

"""
	create_eventfile(logdir, [purge_step=nothing; time=time()]) -> IOStream

Creates a protobuffer events file in the logdir and returns the IO buffer for
writing to it. If `purge_step::Int` is passed then a special event is written
that makes TensorBoard ignore all events before that step (usefull to ignore
part of a calculation that crashed).

Optional keyword argument `prepend` can be passed to prepend a path to the file
name.
"""
function create_eventfile(logdir, purge_step=nothing, time=time(); prepend="")
    hostname = gethostname()
    fname    = prepend*"events.out.tfevents.$time.$hostname"
    fpath    = joinpath(logdir, fname)

    mkpath(dirname(fpath))
    file     = open(fpath, "w")

    # Create the initial log
    if purge_step != nothing
        ev_0 = Event(wall_time=time, step=purge_step, file_version="brain.Event:2")
        write_event(file, ev_0)
        sess_log = TensorBoardLogger.SessionLog(status=TensorBoardLogger.SessionLog_SessionStatus.START)
        ev_0 = Event(wall_time=time, step=purge_step, session_log=sess_log)
        write_event(file, ev_0)
    else
        ev_0 = Event(wall_time=time, step=0, file_version="brain.Event:2")
        write_event(file, ev_0)
    end
	return fname, file
end

"""
	add_eventfile(lg::TBLogger, path::String)

Adds an event file to `lg` with `path` prepended to it's name. It can be used
to create sub-event collection in a single event collection.
"""
function add_eventfile(lg::TBLogger, path="")
    fname, file = create_eventfile(logdir(lg), prepend=path)
    lg.all_files[fname] = file
    return fname
end

# Accessors
"""
    logdir(lg::TBLogger) -> String

Returns the directory where `lg` is writing events data.
"""
logdir(lg::TBLogger)   = lg.logdir

"""
    get_file(lg::TBLogger) -> IOS

Returns the main `file` IOStream object of Logger `lg`.
"""
get_file(lg::TBLogger) = lg.file

"""
    get_file(lg, tags::String...) -> IOS

Returns the `file` IOStream object of Logger `lg` writing to the tag
`tags1/tags2.../tagsN`.
"""
function get_file(lg::TBLogger, tags::String...)
    key = joinpath(tags...)
    if key ∈ keys(lg.all_files)
        return lg.all_files[key]
    else
        return add_eventfile(lg, key)
    end
end

"""
    set_step!(lg, step) -> Int

Sets the iteration counter in the logger to `step`. This counter is used by the
logger when no value is passed by the user.
"""
set_step!(lg::TBLogger, step) = lg.global_step = step

"""
    set_step_increment!(lg, increment) -> Int

Sets the default increment applyed to logger `lg`'s iteration counter
each time logging is performed.

Can be overidden by passing `log_step_increment=some_increment` when logging.
"""
set_step_increment!(lg::TBLogger, Δstep) = lg.step_increment = Δstep

"""
	increment_step!(lg, Δ_Step) -> Int

Increments the step counter in the logger by `Δ_Step` and returns the new value.
"""
increment_step!(lg::TBLogger, Δ_Step) = lg.global_step += Δ_Step

"""
    step(lg) -> Int

Returns the internal step counter of the logger.
"""
step(lg::TBLogger) = lg.global_step

"""
    reset!(lg)

Reset the TBLogger `lg`, deleting everything in it's log directory.
"""
function reset!(lg::TBLogger)
    # close open streams
    for k=keys(lg.all_files)
        close(lg.all_files[k])
        delete!(lg.all_files, k)
    end

    # Overwrite the logdirectoy and create a new base event file
    init_logdir(logdir(lg), tb_overwrite)
    fname   = add_eventfile(lg, "")
    lg.file = get_file(lg, fname)
    set_step!(lg, 0)

    return lg
end

# Implement the AbstractLogger Interface

CoreLogging.catch_exceptions(lg::TBLogger) = false

CoreLogging.min_enabled_level(lg::TBLogger) = lg.min_level

# For now, log everything that is above the lg.min_level
CoreLogging.shouldlog(lg::TBLogger, level, _module, group, id) = true

function CoreLogging.handle_message(lg::TBLogger, level, message, _module, group,
									id, file, line; kwargs...)
    # Unpack the message
    summ    = SummaryCollection()
    i_step = lg.step_increment # :log_step_increment default value

    if !isempty(kwargs)
        data = Vector{Pair{String,Any}}()

        # ∀ (k-v) pairs, decompose values into objects that can be serialized
        for (key,val) in pairs(kwargs)
            # special key describing step increment
            if key == :log_step_increment
                i_step = val
                continue
            end

            preprocess(message*"/$key", val, data)
        end

        # Serialize every object
        for (name,val) in data
            push!(summ.value, summary_impl(name, val))
        end
    end
    iter = increment_step!(lg, i_step)
    write_event(lg.file, make_event(lg, summ, step=iter))
end

######################### Methods for pretty printing ##########################
Base.show(io::IO, tbl::TBLogger) = begin
	str  = "TBLogger(\"$(tbl.logdir)\"), min_level=$(tbl.min_level), "*
		   "purge_step=$(tbl.global_step))"
    Base.print(io, str)
end

Base.show(io::IO, mime::MIME"text/plain", tbl::TBLogger) = begin
	str = """
	TBLogger:
		- Log level     : $(tbl.min_level)
		- Current step  : $(tbl.global_step)
		- Output        : $(tbl.logdir)
		- open files    : $(length(tbl.all_files))
	"""
    Base.print(io, str)
end


text:
"""
    render(x::MyType)

Generic function that defines how a Julia object is rendered.
Should return a `Node` object.

# Examples
```julia
struct MyPlot
    s::Scope
end

WebIO.render(p::MyPlot) = WebIO.render(p.s)
```
"""
function render end
render(x::Union{Node, String}) = x
render(x::Text) = dom"pre"(x.content)
render(::Nothing) = ""
render(x::Any) =
    dom"div"(; setInnerHtml=richest_html(x))

"""
Called after a provider is setup
"""
function re_register_renderables()
    for provider in providers_initialised
        for T in renderable_types
            register_renderable(T, Val(provider))
        end
    end
end

"""
Generic Conversion to Nodes
"""
const mime_order = map(MIME, [ "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/markdown", "application/javascript", "text/plain" ])

function richest_mime(val)
    for mimetype in mime_order
        showable(mimetype, val) && return mimetype
    end
    error("value not writable for any mimetypes")
end

function richest_html(val)
    topmime = richest_mime(val)
    str_repr = stringmime(topmime, val)
    if topmime == MIME("text/html")
        return str_repr # |> WebIO.encode_scripts
    elseif topmime in map(MIME, ["image/png", "image/jpeg"])
        return ""
    elseif topmime == MIME("image/svg+xml")
        return str_repr
    else
        return "
$str_repr
"
    end
end

richest_html(::Nothing) = ""

function htmlstring(val)
    sprint() do io
        show(io, MIME"text/html"(), val)
    end
end

htmlstring(val::AbstractString) = val

"""
Wrap an observable in a scope to enable "live updating."

This method also contains distinct code paths for the cases where the observable
contains "non-simple" data types (in particular, observables that contain
Nodes, Scopes, or Widgets need specially handling).
"""
function observable_to_scope(obs::Observable)
    # Create scope that will contain our output observable for rendering.
    scope = Scope()

    # Create output observable (stored rich html representation of input obs).
    # We must create output via the Observable(scope, ...) constructor so that
    # the onjs call below works (the js is attached to the scope, not the
    # observable itself).

    # Do we need separate code paths for Scope/Widget and Node?
    if isa(obs[], Scope) || isa(obs[], AbstractWidget)
        output = Observable(scope, "obs-scope", render(obs[]))
        map!(render, output, obs)
        ensure_sync(scope, "obs-scope")
        scope.dom = node(ObservableNode(output.id, "obs-scope"))
        return scope
    end

    if isa(obs[], Node)
        output = Observable(
            scope,
            "obs-node",
            obs[],
        )
        map!(output, obs) do value
            if !isa(value, Node)
                @warn "A rendered observable (scopeid=$(scopeid(scope)), obs.id=$(obs.id)) changed from a WebIO node to $(typeof(value))."
                return nothing
            end
            return value
        end
        ensure_sync(scope, "obs-node")
        scope.dom = node(ObservableNode(output.id, "obs-node"))
        return scope
    end

    output = Observable{AbstractString}(
        scope,
        "obs-output",
        richest_html(obs[])
    )

    # Map new values of obs into the output observable and give it a label
    # within the scope.
    map!(output, obs) do data
        # data -> htmlstring(WebIO.render(data)),
        rich = richest_html(data)
        return rich
    end

    # Avoid not-executing

        """
    )
end

"""
A vector of types that are renderable by WebIO.

This exists because for some providers, we need to create some methods when
certain providers are initialized to allow those providers to display custom
types.
This will be removed in the (hopefully near-term) future as we remove providers
out from WebIO and into the appropriate packages.
"""
const renderable_types = Type[]

"""
    @register_renderable(MyType)
    @register_renderable(MyType) do
        # Render definition
    end

Register a type as renderable by WebIO.
This enables your type to be displayed in the appropriate WebIO frontends
(e.g. Jupyter) without any additional work.

This macro may be called either with just the type that you wish to mark as
renderable or with the body of the [`WebIO.render`](@ref) method using do-block
syntax.

The do-block syntax requires parentheses around the typename.
Additionally, due to inconsistencies in the way macros are resolved, the
do-block syntax must be invoked using `@WebIO.register_renderable`
(**not** `WebIO.@register_renderable`).
If the `@WebIO.register_renderable` syntax looks ugly, it might be preferable
to directly import the macro and use it without qualifying its name.

This macro also defines a method for `Base.show` with the `text/html` MIME so
you should not need to define your own.

# Examples
```julia
struct ScatterPlot
    x::Vector{Float64}
    y::Vector{Float64}
end

# Do-block syntax
# Note that the `@` comes before `WebIO`
@WebIO.register_renderable(ScatterPlot) do plot
    # Construct the scatter plot using DOM primitives...
    return node(...)
end

# Do-block syntax with explicit import
using WebIO: @register_renderable
@register_renderable(ScatterPlot) do plot ... end

# Type name syntax
WebIO.render(plot::ScatterPlot) = node(...)
@WebIO.register_renderable ScatterPlot
```
"""
macro register_renderable(typename)
    return register_renderable_macro_helper(esc(typename))
end

macro register_renderable(f::Expr, typename)
    render_method_expr = quote
        local f = $(esc(f))
        WebIO.render(x::$(esc(typename))) = f(x)
    end

    result_expr = register_renderable_macro_helper(esc(typename))
    push!(result_expr.args, render_method_expr)
    return result_expr
end

function register_renderable_macro_helper(
        typename::Union{Symbol, Expr, Type}
)::Expr
    return quote
        push!(renderable_types, $typename)
        function Base.show(
                io::IO,
                m::Union{MIME"text/html", WEBIO_NODE_MIME},
                x::$typename,
        )
            return Base.show(io, m, WebIO.render(x))
        end
    end
end

"""
    register_renderable(MyType)

This function is deprecated. Please use [`WebIO.@register_renderable`](@ref)
instead.

This function was deprecated because it contained too much *magic* (since
*magic* is firmly within the domain of macros).
In particular, this function resorts to `eval`-ing new method definitions for
the types passed into it which is not what a normal function is supposed to do.
"""
function register_renderable(::Type{T}) where T
    Base.depwarn(
        "`WebIO.register_renderable(Type)` is deprecated; use the "
            * "`@WebIO.register_renderable Type` macro instead.",
        :webio_register_renderable_function,
    )
    @eval $(register_renderable_macro_helper(T))
end

function Base.show(io::IO, m::WEBIO_NODE_MIME, node::Node)
    write(io, JSON.json(node))
end

function Base.show(io::IO, m::MIME"text/html", x::Observable)
    show(io, m, WebIO.render(x))
end

function Base.show(io::IO, m::WEBIO_NODE_MIME, x::Union{Observable, AbstractWidget})
    show(io, m, WebIO.render(x))
end

@deprecate render_internal render


text:
using PixelBlocks
using Test

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