text:
using Flux
using Flux: onecold, onehotbatch
using Flux.Losses: logitcrossentropy
using GraphNeuralNetworks
using MLDatasets: Cora
using Statistics, Random
using CUDA
CUDA.allowscalar(false)

function eval_loss_accuracy(X, y, ids, model, g)
    ŷ = model(g, X)
    l = logitcrossentropy(ŷ[:,ids], y[:,ids])
    acc = mean(onecold(ŷ[:,ids] |> cpu) .== onecold(y[:,ids] |> cpu))
    return (loss = round(l, digits=4), acc = round(acc*100, digits=2))
end

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 5f-3             # learning rate
    epochs = 10         # number of epochs
    seed = 17             # set seed > 0 for reproducibility
    usecuda = false      # if true use cuda (if available)
    nhidden = 64        # dimension of hidden features
end

function train(Layer; verbose=false, kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)

    if args.usecuda && CUDA.functional()
        device = Flux.gpu
        args.seed > 0 && CUDA.seed!(args.seed)
    else
        device = Flux.cpu
    end

    # LOAD DATA
    data = Cora.dataset()
    g = GNNGraph(data.adjacency_list) |> device
    X = data.node_features |> device
    y = onehotbatch(data.node_labels, 1:data.num_classes) |> device
    train_ids = data.train_indices |> device
    val_ids = data.val_indices |> device
    test_ids = data.test_indices |> device
    ytrain = y[:,train_ids]

    nin, nhidden, nout = size(X,1), args.nhidden, data.num_classes

    ## DEFINE MODEL
    model = GNNChain(Layer(nin, nhidden),
                    #  Dropout(0.5),
                     Layer(nhidden, nhidden),
                     Dense(nhidden, nout))  |> device

    ps = Flux.params(model)
    opt = ADAM(args.η)

    ## TRAINING
    function report(epoch)
        train = eval_loss_accuracy(X, y, train_ids, model, g)
        test = eval_loss_accuracy(X, y, test_ids, model, g)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
    end

    verbose && report(0)
    for epoch in 1:args.epochs
        gs = Flux.gradient(ps) do
            ŷ = model(g, X)
            logitcrossentropy(ŷ[:,train_ids], ytrain)
        end
        Flux.Optimise.update!(opt, ps, gs)
        verbose && report(epoch)
    end

    train_res = eval_loss_accuracy(X, y, train_ids, model, g)
    test_res = eval_loss_accuracy(X, y, test_ids, model, g)
    return train_res, test_res
end

for Layer in [
            (nin, nout) -> GCNConv(nin => nout, relu),
            (nin, nout) -> GraphConv(nin => nout, relu, aggr=mean),
            (nin, nout) -> SAGEConv(nin => nout, relu),
            (nin, nout) -> GATConv(nin => nout, relu),
            (nin, nout) -> GINConv(Dense(nin, nout, relu), 0.01, aggr=mean),
            (nin, nout) -> ChebConv(nin => nout, 2),
            (nin, nout) -> ResGatedGraphConv(nin => nout, relu),
            # (nin, nout) -> NNConv(nin => nout),  # needs edge features
            # (nin, nout) -> GatedGraphConv(nout, 2),  # needs nin = nout
            # (nin, nout) -> EdgeConv(Dense(2nin, nout, relu)), # Fits the traning set but does not generalize well
              ]

    # @show Layer(2,2)
    train_res, test_res = train(Layer, verbose=false)
    @test train_res.acc > 95
    @test test_res.acc > 70
end


text:
module C2BasisDefinition

"""
C2 basis on the Torus [0,1]
"""

using ..BasisDefinition, ..DynamicDefinition, ValidatedNumerics
import Base: iterate, length
import ..BasisDefinition: one_vector, integral_covector, is_integral_preserving
import ...ResponseManneville: NormKind, Linf

export C2Basis, dual_val, dual_der, C1, C2

struct C1 <: NormKind end
struct C2 <: NormKind end

"""
Equispaced partition of size n of [0,1]
"""
struct EquispacedPartitionInterval{T} <: AbstractVector{T}
	n::Integer
end

function Base.getindex(p::EquispacedPartitionInterval{T}, i::Int) where {T}
	@boundscheck 1 <= i <= p.n+1 || throw(BoundsError(p, i))
	return convert(T, i-1) / p.n
end

EquispacedPartitionInterval(i::Int) = @error "The real type must be specified explicitly"

Base.size(p::EquispacedPartitionInterval) = (p.n+1,)
Base.IndexStyle(::Type{<:EquispacedPartitionInterval}) = IndexLinear()
Base.issorted(p::EquispacedPartitionInterval) = true

struct C2Basis{T<:AbstractVector} <:Basis
	p::T
	# TODO: check in constructor that p is sorted and starts with 0
end
C2Basis(n::Integer) = C2Basis(EquispacedPartitionInterval{Float64}(n))

"""
Return the size of the C2 basisBase.length(S::AverageZero) = length(S.basis)-1
"""
Base.length(b::C2Basis) = 2*length(b.p)

function ϕ(x::Interval{T}) where {T}
    if x ∩ Interval{T}(-1, 1) == ∅
         return zero(x)
    else
        x₋ = x ∩ Interval{T}(-1, 0)
        val₋ = evalpoly(x₋, (1, 0, 0, 10, 15, 6))
        x₊ = x ∩ Interval{T}(0, 1)
        val₊ = evalpoly(x₊, (1, 0, 0, -10, 15, -6))
        return val₋ ∪ val₊
    end
end

function ϕprime(x::Interval{T}) where {T}#Derivative of ϕ
    if x ∩ Interval{T}(-1, 1) == ∅
         return zero(x)
    else
        x₋ = x ∩ Interval{T}(-1, 0)
        val₋ = evalpoly(x₋, (0, 0, 30, 60, 30))
        x₊ = x ∩ Interval{T}(0, 1)
        val₊ = evalpoly(x₊, (0, 0, -30, 60, -30))
        return val₋ ∪ val₊
    end
end

function ν(x::Interval{T}) where {T}
    if x ∩ Interval{T}(-1, 1) == ∅
         return zero(x)
    else
        x₋ = x ∩ Interval{T}(-1, 0)
        val₋ = evalpoly(x₋, (0, 1, 0, -6, -8, -3))
        x₊ = x ∩ Interval{T}(0, 1)
        val₊ = evalpoly(x₊, (0, 1, 0, -6, 8, -3))
        return val₋ ∪ val₊
    end
end

function νprime(x::Interval{T}) where {T} #Derivative of ν

	if x ∩ Interval{T}(-1, 1) == ∅
         return zero(x)
    else
        x₋ = x ∩ Interval{T}(-1, 0)
        val₋ = evalpoly(x₋, (1, 0, -18, -32, -15))
        x₊ = x ∩ Interval{T}(0, 1)
        val₊ = evalpoly(x₊, (1, 0, -18, 32, -15))
        return val₋ ∪ val₊
    end
end

κ(x) = 6*x*(1-x)

function Base.getindex(B::C2Basis, i::Int)
	n = length(B.p)
	@boundscheck 1 <= i <= 2*n || throw(BoundsError(B, i))
	if i<=n
		#println("val")
		return x->ϕ((n-1)*x-i+1)
	elseif n<i<=2*n
        return x->ν((n-1)*x-i+n+1)/(n-1) #we have to subtract n to i
	end
	#return x->κ(x)
end

function basis_element_with_der(B::C2Basis, i::Int)
    n = length(B.p)
	@boundscheck 1 <= i <= 2*n || throw(BoundsError(B, i))
	if i<=n
		#println("val")
		return x->ϕ((n-1)*x-i+1), x-> (n-1)*ϕprime((n-1)*x-i+1)
	elseif n<i<=2*n
        return x->ν((n-1)*x-i+n+1)/(n-1), x-> νprime((n-1)*x-i+n+1) #we have to subtract n to i
	end
	#return x->κ(x) for i in 1:N+1
end

# Check the absolute values in the formula
dual_val(f::Function, fprime::Function, x, der, derder) = f(x)/der
dual_der(f::Function, fprime::Function, x, der, derder) = fprime(x)/der^2-f(x)*derder/der^3
"""
Return (in an iterator) the pairs (i, (x, |T'(x)|)) where x is a preimage of p[i], which
describe the "dual" L* evaluation(p[i])
"""
function Base.iterate(S::DualComposedWithDynamic{T, Dynamic}, state = (1, 1)) where T<:C2Basis
	i, k = state

	if i == length(S.basis)+1
			return nothing
	end

	# remark that this version supposes that for each i there exists a preimage
	# another more specific version should be implemented for maps with
	# incomplete branches

	x = preim(S.dynamic, k, S.basis.p[i], S.ϵ)
    absT′ = abs(der(S.dynamic, x))
    derder = der_der(S.dynamic, x)

    if i <= n
		ret = (x, (f, fprime) -> dual_val(f, fprime, x, der, derder))
	else
		ret = (x, (f, fprime) -> dual_der(f, fprime, x, der, derder))
	end

	#if k == nbranches(S.dynamic)
	#	return ((i, ret), (i+1, 1))
	#else
	#	return ((i, ret), (i, k+1))
	#end
	#if k == nbranches(S.dynamic)
	#	return ((i, (x, absT′, derder)), (i+1, 1))
	#else
	#	return ((i, (x, absT′, derder)), (i, k+1))
	#end
end

BasisDefinition.is_refinement(Bf::C2Basis, Bc::C2Basis) = Bc.p ⊆ Bf.p

function integral_covector(B::C2Basis)
    n = length(B.p)
	return 1/(n-1)*[@interval 0.5; ones(Interval{Float64}, n-2); @interval 0.5; @interval 1//10; zeros(n-2); @interval -1//10]'
end

function one_vector(B::C2Basis)
	return [ones(length(B.p)); zeros(length(B.p))]
end

"""
Return the range of indices of the elements of the basis whose support intersects
with the given dual element (i.e., a pair (y, absT', derder)).
"""
function BasisDefinition.nonzero_on(B::C2Basis, dual_element)
	y = dual_element[1]
    # Note that this cannot rely on arithmetic unless it is verified
    # searchsortedfirst(a, x) return the index of the first value in a greater than or equal to x
	lo = max(searchsortedfirst(B.p, y.lo) -1, 1)
	hi = searchsortedfirst(B.p, y.hi)
	if lo == 1 # 1:N+1 does not make sense and becomes 1:N
		hi = min(hi, length(B))
	end
	return (lo, hi)
end

"""
Given a preimage ```y``` of a point ```x```, this iterator returns
```\\phi_j(y)/T'(y) ```
"""
function Base.iterate(S::ProjectDualElement{T}, state = (S.j_min, :val) ) where T <: C2Basis
    dual = S.dual_element[2]
	j = state[1]
    n = length(S.basis.p)

    if state[2] == :val              # we first fill up the coordinates of the ϕ
		f, fprime = basis_element_with_der(S.basis, j)
		if state[1] == S.j_max
            return ((j, dual(f, fprime)), (S.j_min, :der))
        end
		return ((j, dual(f, fprime)), (j+1, :val))
	end
    if state[2] == :der             # we then fill up the coordinates of the ν
        if state[1] == S.j_max+1
		    return nothing
        end
        f, fprime = basis_element_with_der(S.basis, j+n)
        return ((j+n, dual(f, fprime)), (j+1, :der))
    end
end

using IntervalOptimisation
function infnormoffunction(B::C2Basis, v)
	n = length(B.p)
	maximum = -Inf
	for i in 1:length(B.p)-1
		coeff= v[i]*[1, 0, 0, -10, 15, -6] #coeff for ϕ
    	coeff+= (v[i+n]/(n-1))*[0, 1, 0, -6, 8, -3]  #coeff for ν
    	# coefficients from the right endpoint
    	coeff+= v[i+1]*[0, 0, 0, 10, -15, 6]
		coeff+= (v[i+n+1]/(n-1))*[0, 0, 0, -4, +7, -3]

		dom = Interval(0, 1)
		f(x) = abs(evalpoly(x, coeff))
		maximum =  max(maximum, maximise(f, dom)[1])
	end
	return maximum
end

function infnormofderivative(B::C2Basis, v)
	n = length(B.p)
	maximum = -Inf
	for i in 1:length(B.p)-1
		coeff= (n-1)*v[i]*[0, 0, -30, 60, -30] #coeff for ϕ
    	coeff+= v[i+n]*[1, 0, -18, 32, -15]  #coeff for ν
    	# coefficients from the right endpoint
    	coeff+= (n-1)*v[i+1]*[0, 0, 30, -60, 30]
		coeff+= v[i+n+1]*[0, 0, -12, +28, -15]

		dom = Interval(0, 1)
		f(x) = abs(evalpoly(x, coeff))
		maximum =  max(maximum, maximise(f, dom)[1])
	end
	return maximum
end

C1Norm(B::C2Basis, v) = infnormoffunction(B,v)+infnormofderivative(B,v)

Base.length(S::AverageZero{T}) where T<:C2Basis = length(S.basis)-1

function Base.iterate(S::AverageZero{T}, state = 1) where T<:C2Basis
	n = length(S.basis)÷2
	v = zeros(2*n)

	i = state
	if i == 2*n
		return nothing
	elseif 1 <= i < n
		v[i+1] = 1
		v[1] = -2
	elseif i == n
		v[i+1] = 1
		v[1] = -1/5
	elseif n<i<2*n-1
		v[i+1] = 1
	elseif i==2*n-1
		v[i+1] = 1
		v[1] = 1/5
	end
	return v, state+1
end

end

using RecipesBase

@userplot PlotC2
@recipe function f(h::PlotC2)
	if length(h.args)!= 2 || (typeof(h.args[1])!= C2) || !(typeof(h.args[2])<:AbstractVector)
		error("Plot C2 needs as an input a C2 Basis and a vector")
	end

	B = h.args[1]
	w = h.args[2]

	seriestype := :path
	collect(B), mid.(w)
end

ext: jl

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
#SPG Exit Condition Dict Constructor
import Base.print

export spgExitCondition, print

mutable struct spgExitCondition
    triggered::Nullable{Integer}
    conditions::Dict{Integer,String}
    info::Array{String,1}
end

"""
spgExitCondition(;trigger::Nullable{Integer} = Nullable{Integer}())\n

Use print(ec::spgExitCondition) to see error code and info.\n
"""
function spgExitCondition(;trigger::Nullable{Integer} = Nullable{Integer}())

    vals = ["EXIT_ROOT_FOUND",
            "EXIT_BPSOL_FOUND",
            "EXIT_LEAST_SQUARES",
            "EXIT_OPTIMAL",
            "EXIT_ITERATIONS",
            "EXIT_LINE_ERROR",
            "EXIT_SUBOPTIMAL_BP",
            "EXIT_MATVEC_LIMIT",
            "EXIT_ACTIVE_SET",
            "EXIT_AT_PARETO"]

    info = ["EXIT -- Found a root",
            "EXIT -- Found a BP solution",
            "EXIT -- Found a least-squares solution",
            "EXIT -- Optimal solution found",
            "ERROR EXIT -- Too many iterations",
            "ERROR EXIT -- Linesearch error", #DEVNOTE# Show Line
            "EXIT -- Found a suboptimal BP solution",
            "EXIT -- Maximum matrix-vector operations reached",
            "EXIT -- Found a possible active set",
            "EXIT -- Reached the pareto curve"]

    # Create conditions Dict
    conds = Dict(i => vals[i] for i=1:length(vals))

    # Construct
    exitcon = spgExitCondition(trigger, conds, info)

    return exitcon
end

#DEVNOTE# Not sure if this should be added to base, right now it requires GenSPGL.print()
"""
--------------------------------------------------------------------------------
GenSPGL\n
Use: print(ec::spgExitCondition)\n

Prints the triggered exit condition reached in spgl1.\n
--------------------------------------------------------------------------------
"""
function print(ec::spgExitCondition)

    if isnull(ec.triggered)
        println("Null Trigger. No exit condition reached in spgl1.")
    else
        val = ec.triggered.value
        reason = ec.conditions[val]
        info = ec.info[val]

        # Output
        println("""
        Exit Condition Number:          $(val)\n
        Exit Condition Triggered:       $(reason)\n
        Additional Information:         $(info)

        """)
    end

end


text:
function get_qmf(m::Int)::Tuple{Vector{Float64},Vector{Float64}}
    if m > 10
        throw("unimplemented")
    end

    H,G = Dict{Int,Vector{Float64}}(),Dict{Int,Vector{Float64}}()

    H[1] = Vector{Float64}([0.7071067811865476, 0.7071067811865476])
    H[2] = Vector{Float64}([-0.12940952255126037, 0.2241438680420134, 0.8365163037378079, 0.48296291314453416])
    H[3] = Vector{Float64}([0.03522629188570953, -0.08544127388202666, -0.13501102001025458, 0.45987750211849154, 0.8068915093110925, 0.33267055295008263])
    H[4] = Vector{Float64}([-0.010597401785069032, 0.0328830116668852, 0.030841381835560764, -0.18703481171909309, -0.027983769416859854, 0.6308807679298589, 0.7148465705529157, 0.2303778133088965])
    H[5] = Vector{Float64}([0.0033357252854737712, -0.012580751999081999, -0.006241490212798274, 0.07757149384004572, -0.032244869584638375, -0.24229488706638203, 0.13842814590132074, 0.7243085284377729, 0.6038292697971896, 0.16010239797419293])
    H[6] = Vector{Float64}([-0.0010773010853084796, 0.004777257510945511, 0.0005538422011614961, -0.03158203931748603, 0.027522865530305727, 0.09750160558732304, -0.12976686756726194, -0.22626469396543983, 0.31525035170919763, 0.7511339080210954, 0.49462389039845306, 0.11154074335010947])
    H[7] = Vector{Float64}([0.00035371379997452024, -0.0018016407040474908, 0.0004295779729213665, 0.01255099855609984, -0.01657454163066688, -0.03802993693501441, 0.08061260915108308, 0.07130921926683026, -0.22403618499387498, -0.14390600392856498, 0.4697822874051931, 0.7291320908462351, 0.3965393194819173, 0.07785205408500918])
    H[8] = Vector{Float64}([-0.00011747678412476953, 0.0006754494064505693, -0.00039174037337694705, -0.004870352993451574, 0.008746094047405777, 0.013981027917398282, -0.044088253930794755, -0.017369301001807547, 0.12874742662047847, 0.0004724845739132828, -0.2840155429615469, -0.015829105256349306, 0.5853546836542067, 0.6756307362972898, 0.31287159091429995, 0.05441584224310401])
    H[9] = Vector{Float64}([3.93473203162716e-05, -0.0002519631889427101, 0.00023038576352319597, 0.0018476468830562265, -0.00428150368246343, -0.004723204757751397, 0.022361662123679096, 0.00025094711483145197, -0.06763282906132997, 0.03072568147933338, 0.14854074933810638, -0.09684078322297646, -0.2932737832791749, 0.13319738582500756, 0.6572880780513005, 0.6048231236901112, 0.24383467461259034, 0.038077947363878345])
    H[10] = Vector{Float64}([-1.3264202894521244e-05, 9.358867032006959e-05, -0.00011646685512928545, -0.0006858566949597116, 0.001992405295185056, 0.001395351747052901, -0.010733175483330575, 0.0036065535669561697, 0.033212674059341, -0.029457536821875813, -0.07139414716639708, 0.09305736460357235, 0.12736934033579325, -0.19594627437737705, -0.24984642432731538, 0.2811723436605775, 0.6884590394536035, 0.5272011889317256, 0.1881768000776915, 0.026670057900555554])

    G[1] = Vector{Float64}([-0.7071067811865476, 0.7071067811865476])
    G[2] = Vector{Float64}([-0.48296291314453416, 0.8365163037378079, -0.2241438680420134, -0.12940952255126037])
    G[3] = Vector{Float64}([-0.33267055295008263, 0.8068915093110925, -0.45987750211849154, -0.13501102001025458, 0.08544127388202666, 0.03522629188570953])
    G[4] = Vector{Float64}([-0.2303778133088965, 0.7148465705529157, -0.6308807679298589, -0.027983769416859854, 0.18703481171909309, 0.030841381835560764, -0.0328830116668852, -0.010597401785069032])
    G[5] = Vector{Float64}([-0.16010239797419293, 0.6038292697971896, -0.7243085284377729, 0.13842814590132074, 0.24229488706638203, -0.032244869584638375, -0.07757149384004572, -0.006241490212798274, 0.012580751999081999, 0.0033357252854737712])
    G[6] = Vector{Float64}([-0.11154074335010947, 0.49462389039845306, -0.7511339080210954, 0.31525035170919763, 0.22626469396543983, -0.12976686756726194, -0.09750160558732304, 0.027522865530305727, 0.03158203931748603, 0.0005538422011614961, -0.004777257510945511, -0.0010773010853084796])
    G[7] = Vector{Float64}([-0.07785205408500918, 0.3965393194819173, -0.7291320908462351, 0.4697822874051931, 0.14390600392856498, -0.22403618499387498, -0.07130921926683026, 0.08061260915108308, 0.03802993693501441, -0.01657454163066688, -0.01255099855609984, 0.0004295779729213665, 0.0018016407040474908, 0.00035371379997452024])
    G[8] = Vector{Float64}([-0.05441584224310401, 0.31287159091429995, -0.6756307362972898, 0.5853546836542067, 0.015829105256349306, -0.2840155429615469, -0.0004724845739132828, 0.12874742662047847, 0.017369301001807547, -0.044088253930794755, -0.013981027917398282, 0.008746094047405777, 0.004870352993451574, -0.00039174037337694705, -0.0006754494064505693, -0.00011747678412476953])
    G[9] = Vector{Float64}([-0.038077947363878345, 0.24383467461259034, -0.6048231236901112, 0.6572880780513005, -0.13319738582500756, -0.2932737832791749, 0.09684078322297646, 0.14854074933810638, -0.03072568147933338, -0.06763282906132997, -0.00025094711483145197, 0.022361662123679096, 0.004723204757751397, -0.00428150368246343, -0.0018476468830562265, 0.00023038576352319597, 0.0002519631889427101, 3.93473203162716e-05])
    G[10] = Vector{Float64}([-0.026670057900555554, 0.1881768000776915, -0.5272011889317256, 0.6884590394536035, -0.2811723436605775, -0.24984642432731538, 0.19594627437737705, 0.12736934033579325, -0.09305736460357235, -0.07139414716639708, 0.029457536821875813, 0.033212674059341, -0.0036065535669561697, -0.010733175483330575, -0.001395351747052901, 0.001992405295185056, 0.0006858566949597116, -0.00011646685512928545, -9.358867032006959e-05, -1.3264202894521244e-05])

    return H[m],G[m]
end

"""
Discrete wavelet transform of vector

:param n: size of x
:type n: Int
:param x: vector to be transformed
:type x: AbstractVector{Float64}
:param m: number of vanishing moments
:type m: Int
:param H: low pass filter
:type H: Vector{Float64}
:param G: high pass filter
:type G: Vector{Float64}
"""
function dwtvec(n::Int,
    x::AbstractVector{Float64},
    m::Int,
    H::Vector{Float64},
    G::Vector{Float64},
    )::Tuple{AbstractVector{Float64},AbstractVector{Float64}}
    d = [sum(G[k]*x[1+mod(2*(i-1)+k-1,n)] for k in 1:2*m) for i in 1:n÷2]
    s = [sum(H[k]*x[1+mod(2*(i-1)+k-1,n)] for k in 1:2*m) for i in 1:n÷2]
    if issparse(x)
        d,s = sparse(d),sparse(s)
    end
    return d,s
end

"""
Inverse wavelet transform of vector

:param n: size of d and s
:type n: Int
:param d: difference vector
:type d: AbstractVector{Float64}
:param s: average vector
:type s: AbstractVector{Float64}
:param m: number of vanishing moments
:type m: Int
:param H: low pass filter
:type H: Vector{Float64}
:param G: high pass filter
:type G: Vector{Float64}
"""
function idwtvec(n::Int,
    d::AbstractVector{Float64},
    s::AbstractVector{Float64},
    m::Int,
    H::Vector{Float64},
    G::Vector{Float64},
    )::AbstractVector{Float64}
    x = zeros(Float64,2*n)
    for i in 1:n
        x[2*(i-1)+1] = sum(H[2*(k-1)+1]*s[1+mod(i-k,n)]+G[2*(k-1)+1]*d[1+mod(i-k,n)] for k in 1:m)
        x[2*(i-1)+2] = sum(H[2*(k-1)+2]*s[1+mod(i-k,n)]+G[2*(k-1)+2]*d[1+mod(i-k,n)] for k in 1:m)
    end
    if issparse(d) & issparse(s)
        x = sparse(x)
    end
    return x
end

"""
Wavelet transform of a matrix

wt(X) -> | A C |
         | B T |

:param n: number of rows/columns in X
:type n: Int
:param X: matrix to be transformed
:param m: number of vanishing moments
:type m: Int
:param H: low pass filter
:type H: Vector{Float64}
:param G: high pass filter
:type G: Vector{Float64}

:return: A,B,C,T blocks
"""
function dwtmat(n::Int,
    X::AbstractMatrix{Float64},
    m::Int,
    H::Vector{Float64},
    G::Vector{Float64},
    )

    if size(X) != (n,n)
        throw("invalid matrix size")
    end

    if n%2 != 0
        throw("size of matrix not divisible by 2")
    end

    if length(H) != 2*m || length(G) != 2*m
        throw("invalid quadrature mirror filters")
    end

    if issparse(X)
        Y = spzeros(size(X,1),size(X,2))
    else
        Y = Matrix{Float64}(undef,size(X))
    end

    for i in 1:n
        Y[i,1:n÷2],Y[i,n÷2+1:n] = dwtvec(n,X[i,:],m,H,G)
    end

    for i in 1:n
        Y[1:n÷2,i],Y[n÷2+1:n,i] = dwtvec(n,Y[:,i],m,H,G)
    end

    return Y[1:n÷2,1:n÷2],Y[n÷2+1:n,1:n÷2],Y[1:n÷2,n÷2+1:n],Y[n÷2+1:n,n÷2+1:n]
end

"""
Inverse wavelet transform of a matrix

:param n: number of rows/columns in X
:type n: Int
:param A: A block
:type A:: AbstractMatrix{Float64}
:param B: B block
:type B:: AbstractMatrix{Float64}
:param C: C block
:type C:: AbstractMatrix{Float64}
:param T: T block
:type T:: AbstractMatrix{Float64}
:param m: number of vanishing moments
:type m: Int
:param H: low pass filter
:type H: Vector{Float64}
:param G: high pass filter
:type G: Vector{Float64}

:return: A,B,C,T blocks
"""
function idwtmat(n::Int,
    m::Int,
    H::Vector{Float64},
    G::Vector{Float64},
    A::AbstractMatrix{Float64},
    B::AbstractMatrix{Float64},
    C::AbstractMatrix{Float64},
    T::AbstractMatrix{Float64},
    )::AbstractMatrix{Float64}

    if issparse(A) & issparse(B) & issparse(C) & issparse(T)
        Y = spzeros(Float64,n*2,n*2)
    else
        Y = Matrix{Float64}(undef,n*2,n*2)
    end

    for i in 1:n
        Y[:,i] = idwtvec(n,A[:,i],B[:,i],m,H,G)
        Y[:,n+i] = idwtvec(n,C[:,i],T[:,i],m,H,G)
    end

    for i in 1:2*n
        Y[i,:] = idwtvec(n,Y[i,1:n],Y[i,n+1:2*n],m,H,G)
    end

    return Y
end

"""
Wavelet transform of a matrix of size (n^2, n^2). The matrix is obatined
by discretizing a 2 dimensional differential operator on a n-by-n grid.

:param n: size of the grid
:type n: Int
:param X: matrix of size (n^2, n^2)
:type X:: AbstractMatrix{Float64}
:param m: number of vanishing moments
:type m: Int
:param H: low pass filter
:type H: Vector{Float64}
:param G: high pass filter
:type G: Vector{Float64}

:return: wavelet transform
"""
function dwtmat2d(n::Int,
    X::AbstractMatrix{Float64},
    m::Int,
    H::Vector{Float64},
    G::Vector{Float64},
    )
    # Check inputs
    if size(X) != (n^2,n^2)
        throw("invalid matrix size of $(size(X))")
    end

    if n%2 != 0
        throw("size of matrix not divisible by 2")
    end

    if length(H) != 2*m || length(G) != 2*m
        throw("invalid quadrature mirror filters")
    end

    Y = copy(X)
    k = n÷2

    # Transform rows of n consecutive entries
    for i in 1:n^2
        for j in 1:n
            d,s = dwtvec(n,Y[i,(j-1)*n+1:j*n],m,H,G)
            Y[i,(j-1)*n+1:(j-1)*n+k] = d
            Y[i,(j-1)*n+k+1:j*n] = s
        end
    end

    # Transform columns of n consecutive entries
    for i in 1:n^2
        for j in 1:n
            d,s = dwtvec(n,Y[(j-1)*n+1:j*n,i],m,H,G)
            Y[(j-1)*n+1:(j-1)*n+k,i] = d
            Y[(j-1)*n+k+1:j*n,i] = s
        end
    end

    # Transform rows of n entries with a stride of n
    for i in 1:n^2
        for j in 1:n
            d,s = dwtvec(n,Y[i,j:n:(n-1)*n+j],m,H,G)
            Y[i,j:n:(k-1)*n+j] = d
            Y[i,k*n+j:n:(n-1)*n+j] = s
        end
    end

    # Transform column of n entries with a stride of n
    for i in 1:n^2
        for j in 1:n
            d,s = dwtvec(n,Y[j:n:(n-1)*n+j,i],m,H,G)
            Y[j:n:(k-1)*n+j,i] = d
            Y[k*n+j:n:(n-1)*n+j,i] = s
        end
    end

    # Re-indexing
    grid = reshape(collect(1:n^2),(n,n))
    index = vcat(collect(Iterators.flatten(grid[1:n÷2,1:n÷2])),
                collect(Iterators.flatten(grid[n÷2+1:n,1:n÷2])),
                collect(Iterators.flatten(grid[1:n÷2,n÷2+1:n])),
                collect(Iterators.flatten(grid[n÷2+1:n,n÷2+1:n])),
            )

    return Y[index[1:3*k^2],index[1:3*k^2]],Y[index[3*k^2+1:n^2],index[1:3*k^2]],
        Y[index[1:3*k^2],index[3*k^2+1:n^2]],Y[index[3*k^2+1:n^2],index[3*k^2+1:n^2]]
end

function solve(A::Matrix{Float64},
    b::Matrix{Float64},
    )::Matrix{Float64}

    if size(A) != (size(b)[1]*size(b)[2],size(b)[1]*size(b)[2])
        throw(DimensionMismatch("size of A not equal to size of b"))
    end

    n = size(A)[1]
    F = LinearAlgebra.svd(A)
    sig = Vector{Float64}(undef,n)
    for i in 1:n
        if F.S[i] > 10e-7
            sig[i] = 1.0/F.S[i]
        else
            sig[i] = 0.0
        end
    end
    println("coarsest scale condition number $(F.S[1]/F.S[n-1]) and $(F.S[n])")

    x = reshape(transpose(F.Vt)*(sig.*(transpose(F.U)*reshape(b,n,1))),size(b))

    return x
end

max_stars_repo_name: xinshuoyang/FactNSF

max_issues_repo_name: xinshuoyang/FactNSF

max_forks_repo_name: xinshuoyang/FactNSF


text:
"""
$(TYPEDEF)
"""
abstract type BasicEnsembleAlgorithm <: EnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleThreads <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleDistributed <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleSplitThreads <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleSerial <: BasicEnsembleAlgorithm end

#=
if (kwargs[:parallel_type] == != :none && kwargs[:parallel_type] == != :threads)
  error("Distributed arrays cannot be generated via none or threads")
end
(batch_size != trajectories) && warn("batch_size and reductions are ignored when !collect_result")

elapsed_time = @elapsed u = DArray((trajectories,)) do II
    solve_batch(prob,alg,kwargs[:parallel_type] ==,II[1],pmap_batch_size,kwargs...)
end
return EnsembleSolution(u,elapsed_time,false)
=#

@noinline function parallel_type_warn()
  @warn "parallel_type has been deprecated. Please refer to the docs for the new dispatch-based system."
end

@noinline function not_recognized_error()
  @error "parallel_type value not recognized"
end

@noinline function num_monte_warn()
  @warn "num_monte has been replaced by trajectories"
end

function __solve(prob::AbstractEnsembleProblem,
                 alg::Union{DEAlgorithm,Nothing};
                 kwargs...)
    #=
    if alg isa EnsembleAlgorithm
      @error "You forgot to pass a DE solver algorithm! Only a EnsembleAlgorithm has been supplied. Exiting"
    end
    =#
    if :parallel_type ∈ keys(kwargs)
      if kwargs[:parallel_type] == :none
        ensemblealg = EnsembleSerial()
      elseif kwargs[:parallel_type] == :pmap || kwargs[:parallel_type] == :parfor
        ensemblealg = EnsembleDistributed()
      elseif kwargs[:parallel_type] == :threads
        ensemblealg = EnsembleThreads()
      elseif kwargs[:parallel_type] == :split_threads
        ensemblealg = EnsembleSplitThreads()
      else
        not_recognized_error()
      end
    elseif alg isa EnsembleAlgorithm
      # Assume DifferentialEquations.jl is being used, so default alg
      ensemblealg = alg
      alg = nothing
    else
      ensemblealg = EnsembleThreads()
    end
    if :num_monte ∈ keys(kwargs)
      num_monte_warn()
      trajectories = kwargs[:num_monte]
    else
      @assert :trajectories ∈ keys(kwargs)
      trajectories = kwargs[:trajectories]
    end
    if :parallel_type ∈ keys(kwargs)
      parallel_type_warn()
    end
    __solve(prob,alg,ensemblealg;trajectories=trajectories,kwargs...)
end

@noinline function rerun_warn()
  @warn("output_func should return (out,rerun). See docs for updated details")
end
tighten_container_eltype(u::Vector{Any}) = map(identity, u)
tighten_container_eltype(u) = u

function __solve(prob::AbstractEnsembleProblem,
                 alg::Union{DEAlgorithm,Nothing},
                 ensemblealg::BasicEnsembleAlgorithm;
                 trajectories, batch_size = trajectories,
                 pmap_batch_size = batch_size÷100 > 0 ? batch_size÷100 : 1, kwargs...)

  num_batches = trajectories ÷ batch_size
  num_batches < 1 && error("trajectories ÷ batch_size cannot be less than 1, got $num_batches")
  num_batches * batch_size != trajectories && (num_batches += 1)

  if num_batches == 1 && prob.reduction === DEFAULT_REDUCTION
    elapsed_time = @elapsed u = solve_batch(prob,alg,ensemblealg,1:trajectories,pmap_batch_size;kwargs...)
    _u = tighten_container_eltype(u)
    return EnsembleSolution(_u,elapsed_time,true)
  end

  converged::Bool = false
  i = 1
  II = (batch_size*(i-1)+1):batch_size*i

  batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...)

  u = prob.u_init === nothing ? similar(batch_data, 0) : prob.u_init
  u,converged = prob.reduction(u,batch_data,II)
  elapsed_time = @elapsed for i in 2:num_batches
    converged && break
    if i == num_batches
      II = (batch_size*(i-1)+1):trajectories
    else
      II = (batch_size*(i-1)+1):batch_size*i
    end
    batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...)
    u,converged = prob.reduction(u,batch_data,II)
  end

  _u = tighten_container_eltype(u)

  return EnsembleSolution(_u,elapsed_time,converged)

end

function batch_func(i,prob,alg;kwargs...)
  iter = 1
  _prob = prob.safetycopy ? deepcopy(prob.prob) : prob.prob
  new_prob = prob.prob_func(_prob,i,iter)
  rerun = true
  x = prob.output_func(solve(new_prob,alg;kwargs...),i)
  if !(typeof(x) <: Tuple)
      rerun_warn()
      _x = (x,false)
  else
    _x = x
  end
  rerun = _x[2]
  while rerun
      iter += 1
      _prob = prob.safetycopy ? deepcopy(prob.prob) : prob.prob
      new_prob = prob.prob_func(_prob,i,iter)
      x = prob.output_func(solve(new_prob,alg;kwargs...),i)
      if !(typeof(x) <: Tuple)
          rerun_warn()
          _x = (x,false)
      else
        _x = x
      end
      rerun = _x[2]
  end
  _x[1]
end

function solve_batch(prob,alg,ensemblealg::EnsembleDistributed,II,pmap_batch_size;kwargs...)
  wp=CachingPool(workers())

  # Fix the return type of pmap
  #=
  function f()
    batch_func(1,prob,alg;kwargs...)
  end
  T = Core.Compiler.return_type(f,Tuple{})
  =#

  batch_data = pmap(wp,II,batch_size=pmap_batch_size) do i
    batch_func(i,prob,alg;kwargs...)
  end

  tighten_container_eltype(batch_data)
end

function responsible_map(f,II...)
  batch_data = []
  sizehint!(batch_data,length(II[1]))
  for i in 1:length(II[1])
    @inbounds push!(batch_data, f(getindex.(II,i)...))
  end
  identity.(batch_data)
end

function SciMLBase.solve_batch(prob,alg,::EnsembleSerial,II,pmap_batch_size;kwargs...)
  batch_data = responsible_map(II) do i
    SciMLBase.batch_func(i,prob,alg;kwargs...)
  end
  SciMLBase.tighten_container_eltype(batch_data)
end

function solve_batch(prob,alg,ensemblealg::EnsembleThreads,II,pmap_batch_size;kwargs...)

  if length(II) == 1 || Threads.nthreads() == 1
    return solve_batch(prob,alg,EnsembleSerial(),II,pmap_batch_size;kwargs...)
  end

  if typeof(prob.prob) <: AbstractJumpProblem && length(II) != 1
    probs = [deepcopy(prob.prob) for i in 1:Threads.nthreads()]
  else
    probs = prob.prob
  end

  #
  batch_size = length(II)÷Threads.nthreads()

  batch_data = tmap(1:Threads.nthreads()) do i
    if i == Threads.nthreads()
      I_local = II[(batch_size*(i-1)+1):end]
    else
      I_local = II[(batch_size*(i-1)+1):(batch_size*i)]
    end
    solve_batch(prob,alg,EnsembleSerial(),I_local,pmap_batch_size;kwargs...)
  end
  tighten_container_eltype(batch_data)
end

function tmap(f,args...)
  batch_data = Vector{Core.Compiler.return_type(f,Tuple{typeof.(getindex.(args,1))...})}(undef,length(args[1]))
  Threads.@threads for i in 1:length(args[1])
      batch_data[i] = f(getindex.(args,i)...)
  end
  batch_data
end

function solve_batch(prob,alg,::EnsembleSplitThreads,II,pmap_batch_size;kwargs...)
  wp=CachingPool(workers())
  N = nworkers()
  batch_size = length(II)÷N

  # Fix the return type of pmap
  #=
  function f()
    i = 1
    I_local = II[(batch_size*(i-1)+1):end]
    solve_batch(prob,alg,EnsembleThreads(),I_local,pmap_batch_size;kwargs...)
  end
  T = Core.Compiler.return_type(f,Tuple{})
  =#

  batch_data = let
    pmap(wp,1:N,batch_size=pmap_batch_size) do i
      if i == N
        I_local = II[(batch_size*(i-1)+1):end]
      else
        I_local = II[(batch_size*(i-1)+1):(batch_size*i)]
      end
      solve_batch(prob,alg,EnsembleThreads(),I_local,pmap_batch_size;kwargs...)
    end
  end
  reduce(vcat,batch_data)
end


text:
abstract InferenceAlgorithm{P}
abstract Sampler{T<:InferenceAlgorithm}

## Fallback functions

Base.run(spl :: Sampler) = error("[sample]: unmanaged inference algorithm: $(typeof(spl))")

assume(spl, distr :: Distribution) =
  error("[assume]: unmanaged inference algorithm: $(typeof(spl))")

observe(spl, weight :: Float64) =
  error("[observe]: unmanaged inference algorithm: $(typeof(spl))")

predict(spl, var_name :: Symbol, value) =
  error("[predict]: unmanaged inference algorithm: $(typeof(spl))")

## Default functions
function sample(model::Function, alg :: InferenceAlgorithm)
  global sampler = ParticleSampler{typeof(alg)}(model, alg);
  Base.run(sampler)
end

assume(spl :: Sampler, distr :: Distribution)  = rand(current_trace(), distr)
observe(spl :: Sampler, score :: Float64)      = produce(score)

function predict(spl :: Sampler, v_name :: Symbol, value)
  task = current_task()
  if haskey(task.storage, :turing_predicts)
    predicts = task.storage[:turing_predicts]
  else
    predicts = Dict{Symbol,Any}()
  end
  predicts[v_name] = value
  task.storage[:turing_predicts] = predicts
end


text:
module CameraApp

using VideoIO
using ImageInTerminal

export demo

function demo()
    ImageInTerminal.use_24bit()
    f = opencamera()
    outH, outW = displaysize(stdout)
    try
        print("\x1b[?25l") # hide cursor
        while true
            buf = IOBuffer()
            img = read(f)
            outH, outW = displaysize(stdout)
            print(buf, "\x1b[999D\x1b[$(outH)A")
            ImageInTerminal.imshow(buf, img, (outH, outW))
            println(buf |> take! |> String)
            sleep(0.01)
        end
    catch e
        close(f)
        print("\x1b[?25h") # unhide cursor
    finally
        println("exit")
    end
end

function real_main()
    # Write julia logic what you want
    demo()
end

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

end # module


text:
# Julia wrappers and extensions for Hashpipe header fitshead.h
# Automatically generated using Clang.jl before manual code extensions

function hgeti2(hstring, keyword, val)
    ccall((:hgeti2, libhashpipe), Cint, (Cstring, Cstring, Ptr{Int16}), hstring, keyword, val)
end

function hgeti4c(hstring, keyword, wchar, val)
    ccall((:hgeti4c, libhashpipe), Cint, (Cstring, Cstring, Cstring, Ptr{Cint}), hstring, keyword, wchar, val)
end

function hgeti4(hstring, keyword, val)
    ccall((:hgeti4, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cint}), hstring, keyword, val)
end

function hgetu4(hstring, keyword, val)
    ccall((:hgetu4, libhashpipe), Cint, (Cstring, Cstring, Ptr{UInt32}), hstring, keyword, val)
end

function hgetr4(hstring, keyword, val)
    ccall((:hgetr4, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cfloat}), hstring, keyword, val)
end

function hgetr8c(hstring, keyword, wchar, val)
    ccall((:hgetr8c, libhashpipe), Cint, (Cstring, Cstring, Cstring, Ptr{Cdouble}), hstring, keyword, wchar, val)
end

function hgetr8(hstring, keyword, val)
    ccall((:hgetr8, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cdouble}), hstring, keyword, val)
end

function hgetra(hstring, keyword, ra)
    ccall((:hgetra, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cdouble}), hstring, keyword, ra)
end

function hgetdec(hstring, keyword, dec)
    ccall((:hgetdec, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cdouble}), hstring, keyword, dec)
end

function hgetdate(hstring, keyword, date)
    ccall((:hgetdate, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cdouble}), hstring, keyword, date)
end

function hgetl(hstring, keyword, lval)
    ccall((:hgetl, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cint}), hstring, keyword, lval)
end

function hgetsc(hstring, keyword, wchar, lstr, string)
    ccall((:hgetsc, libhashpipe), Cint, (Cstring, Cstring, Cstring, Cint, Cstring), hstring, keyword, wchar, lstr, string)
end

function hgets(hstring, keyword, lstr, string)
    ccall((:hgets, libhashpipe), Cint, (Cstring, Cstring, Cint, Cstring), hstring, keyword, lstr, string)
end

function hgetndec(hstring, keyword, ndec)
    ccall((:hgetndec, libhashpipe), Cint, (Cstring, Cstring, Ptr{Cint}), hstring, keyword, ndec)
end

function hgetc(hstring, keyword)
    ccall((:hgetc, libhashpipe), Cstring, (Cstring, Cstring), hstring, keyword)
end

function hgetc_thread_safe(hstring, keyword, value_buffer)
    ccall((:hgetc_thread_safe, libhashpipe), Cstring, (Cstring, Cstring, Cstring), hstring, keyword, value_buffer)
end

function ksearch(hstring, keyword)
    ccall((:ksearch, libhashpipe), Cstring, (Cstring, Cstring), hstring, keyword)
end

function blsearch(hstring, keyword)
    ccall((:blsearch, libhashpipe), Cstring, (Cstring, Cstring), hstring, keyword)
end

function strsrch(s1, s2)
    ccall((:strsrch, libhashpipe), Cstring, (Cstring, Cstring), s1, s2)
end

function strnsrch(s1, s2, ls1)
    ccall((:strnsrch, libhashpipe), Cstring, (Cstring, Cstring, Cint), s1, s2, ls1)
end

function strcsrch(s1, s2)
    ccall((:strcsrch, libhashpipe), Cstring, (Cstring, Cstring), s1, s2)
end

function strncsrch(s1, s2, ls1)
    ccall((:strncsrch, libhashpipe), Cstring, (Cstring, Cstring, Cint), s1, s2, ls1)
end

function hlength(header, lhead)
    ccall((:hlength, libhashpipe), Cint, (Cstring, Cint), header, lhead)
end

function gethlength(header)
    ccall((:gethlength, libhashpipe), Cint, (Cstring,), header)
end

function str2ra(in)
    ccall((:str2ra, libhashpipe), Cdouble, (Cstring,), in)
end

function str2dec(in)
    ccall((:str2dec, libhashpipe), Cdouble, (Cstring,), in)
end

function isnum(string)
    ccall((:isnum, libhashpipe), Cint, (Cstring,), string)
end

function notnum(string)
    ccall((:notnum, libhashpipe), Cint, (Cstring,), string)
end

function numdec(string)
    ccall((:numdec, libhashpipe), Cint, (Cstring,), string)
end

function strfix(string, fillblank, dropzero)
    ccall((:strfix, libhashpipe), Cvoid, (Cstring, Cint, Cint), string, fillblank, dropzero)
end

function getltime()
    ccall((:getltime, libhashpipe), Cstring, ())
end

function getutime()
    ccall((:getutime, libhashpipe), Cstring, ())
end

function hputi2(hstring, keyword, ival)
    ccall((:hputi2, libhashpipe), Cint, (Cstring, Cstring, Int16), hstring, keyword, ival)
end

function hputi4(hstring, keyword, ival)
    ccall((:hputi4, libhashpipe), Cint, (Cstring, Cstring, Cint), hstring, keyword, ival)
end

# Auto-convert julia string/int to Cstring/Cint
function hputi4(p_hstring::Ptr{UInt8}, p_keyword::String, p_ival::Int)
    error::Int = ccall((:hputi4, libhashpipestatus),
                    Int, (Ptr{UInt8}, Cstring, Cint),
                    p_hstring, Cstring(pointer(p_keyword)), Cint(p_ival))
    return error
end

function hputu4(hstring, keyword, ival)
    ccall((:hputu4, libhashpipe), Cint, (Cstring, Cstring, UInt32), hstring, keyword, ival)
end

function hputr4(hstring, keyword, rval)
    ccall((:hputr4, libhashpipe), Cint, (Cstring, Cstring, Cfloat), hstring, keyword, rval)
end

function hputr8(hstring, keyword, dval)
    ccall((:hputr8, libhashpipe), Cint, (Cstring, Cstring, Cdouble), hstring, keyword, dval)
end

function hputnr8(hstring, keyword, ndec, dval)
    ccall((:hputnr8, libhashpipe), Cint, (Cstring, Cstring, Cint, Cdouble), hstring, keyword, ndec, dval)
end

function hputs(hstring, keyword, cval)
    ccall((:hputs, libhashpipe), Cint, (Cstring, Cstring, Cstring), hstring, keyword, cval)
end

# Auto-convert Julia string to Cstring
function hputs(p_hstring::Ptr{UInt8}, p_keyword::String, p_cval::String)
    error::Int = ccall((:hputs, libhashpipestatus),
                    Int, (Ptr{UInt8}, Cstring, Cstring),
                    p_hstring, Cstring(pointer(p_keyword)), Cstring(pointer(p_cval)))
    return error
end

function hputm(hstring, keyword, cval)
    ccall((:hputm, libhashpipe), Cint, (Cstring, Cstring, Cstring), hstring, keyword, cval)
end

function hputcom(hstring, keyword, comment)
    ccall((:hputcom, libhashpipe), Cint, (Cstring, Cstring, Cstring), hstring, keyword, comment)
end

function hputra(hstring, keyword, ra)
    ccall((:hputra, libhashpipe), Cint, (Cstring, Cstring, Cdouble), hstring, keyword, ra)
end

function hputdec(hstring, keyword, dec)
    ccall((:hputdec, libhashpipe), Cint, (Cstring, Cstring, Cdouble), hstring, keyword, dec)
end

function hputl(hstring, keyword, lval)
    ccall((:hputl, libhashpipe), Cint, (Cstring, Cstring, Cint), hstring, keyword, lval)
end

function hputc(hstring, keyword, cval)
    ccall((:hputc, libhashpipe), Cint, (Cstring, Cstring, Cstring), hstring, keyword, cval)
end

function hdel(hstring, keyword)
    ccall((:hdel, libhashpipe), Cint, (Cstring, Cstring), hstring, keyword)
end

function hadd(hplace, keyword)
    ccall((:hadd, libhashpipe), Cint, (Cstring, Cstring), hplace, keyword)
end

function hchange(hstring, keyword1, keyword2)
    ccall((:hchange, libhashpipe), Cint, (Cstring, Cstring, Cstring), hstring, keyword1, keyword2)
end

function ra2str(string, lstr, ra, ndec)
    ccall((:ra2str, libhashpipe), Cvoid, (Cstring, Cint, Cdouble, Cint), string, lstr, ra, ndec)
end

function dec2str(string, lstr, dec, ndec)
    ccall((:dec2str, libhashpipe), Cvoid, (Cstring, Cint, Cdouble, Cint), string, lstr, dec, ndec)
end

function deg2str(string, lstr, deg, ndec)
    ccall((:deg2str, libhashpipe), Cvoid, (Cstring, Cint, Cdouble, Cint), string, lstr, deg, ndec)
end

function num2str(string, num, field, ndec)
    ccall((:num2str, libhashpipe), Cvoid, (Cstring, Cdouble, Cint, Cint), string, num, field, ndec)
end

function setheadshrink(hsh)
    ccall((:setheadshrink, libhashpipe), Cvoid, (Cint,), hsh)
end

function setleaveblank(hsh)
    ccall((:setleaveblank, libhashpipe), Cvoid, (Cint,), hsh)
end


text:
using LinearAlgebra

using Flux
using Flux.Zygote: dropgrad

import Statistics: mean

"""
    smooth_l1_loss(y, fx; δ)

l1 loss w/ clipping. Also known as the huber loss. Only usable w/ numbers (see `huber_loss` for vector form)
"""
function smooth_l1_loss(y::Number, fx::Number, δ=1)
    α = abs(y - fx)
    α <= δ && return 0.5f0 * α ^ 2
    δ * α - (0.5f0 * δ ^ 2)
end

"""
    huber_loss(y, fx; δ)

Huber loss. Convenience for smoot_l1_loss
"""
function huber_loss(y, fx; δ=1)
    return smooth_l1_loss.(y, fx, δ)
end

"""
    mean_huber_loss(y, fx; δ)

mean huber loss.
"""
function mean_huber_loss(y, fx; δ=1)
    return mean(huber_loss(y, fx; δ=δ))
end

get_cart_idx(a, l) = CartesianIndex.(a, 1:l)
Flux.Zygote.@nograd get_cart_idx

function loss(lu::QLearning, model, s_t, a_t, s_tp1, r, terminal, target_model)

    q_tp1 = if target_model isa Nothing
        maximum(model(s_tp1); dims=1)[1, :]
    else
        maximum(target_model(s_tp1); dims=1)[1, :]
    end
    targets = dropgrad(r + lu.γ*(one(lu.γ) .- terminal).*q_tp1)

    action_idx = get_cart_idx(a_t, length(terminal))
    q_t = model(s_t)[action_idx]
    return lu.loss(q_t, targets)
end

function loss(lu::DoubleQLearning, model, s_t, a_t, s_tp1, r, terminal, target_model)
    γ = lu.γ.*(1 .- terminal)

    action_idx = get_cart_idx(a_t, length(terminal))

    q̃_tp1 = dropgrad(model(s_tp1))
    q̃_tp1_argmax = dropgrad(findmax(q̃_tp1; dims=1))
    action_tp1 = [q̃_tp1_argmax[2][i] for i in 1:length(terminal)]

    q_tp1 = if target_model isa Nothing
        dropgrad(model(s_tp1)[action_tp1])
    else
        dropgrad(target_model(s_tp1)[action_tp1])
    end

    target = dropgrad(r .+ γ.*q_tp1)
    q_t = model(s_t)[action_idx]

    return lu.loss(q_t, target)
end

struct TDLearning end

# function loss(lu::TDLearning, model, s_t, a_t, s_tp1, r, terminal, target_model, horde::H) where {H<:Horde}

#     # get GVF horde parameters

#     p = [RLCore.get(horde, s_t[:,i], a_t, s_tp1[:,i]) for i in 1:length(terminal)]

#     c = hcat(getindex.(p, 1)...)
#     γ = hcat(getindex.(p, 2)...)
#     π = hcat(getindex.(p, 3)...)

#     v_t = model(s_t)
#     v_tp1 = target_model(s_tp1)

#     target = dropgrad(c .+ γ.*v_tp1) # -> Matrix (preds × batch_size)

#     return sum((target .- v_t).^2) * (1 // length(terminal))
# end

# loss(lu::TDLearning, model, s_t, a_t, s_tp1, r, terminal, target_model::Nothing, horde::H) where {H<:Horde} =
#     loss(lu, model, s_t, a_t, s_tp1, r, terminal, (x)->Flux.data(model(x)), horde)


text:
"""
    A = JopFft(T, n...[; dims=()])

where `A` is a Fourier transform, the domain of A is (T,n).  If `dims`
is not set, then `A` is an N-dimensional Fourier transform and where `N=ndims(dom)`.
Otherwise, `A` is a Fourier transform over the dimensions specifed by `dims`.

For a complex-to-complex transform: `range(A)::JetSpace`.  For a real-to-complex
transform: `range(A)::JetSpaceSymmetric`.  Normally, JopFft is responsible for constructing
range(A).  In the event that you need to manually construct this space, we provide
a convenience method:

    R = symspace(JopLn{JopFft_df!}, T, symdim, n...)

where `T` is the precision (generally `Float32` or `Float64`), `n...` is
the dimensions of the domain of `A`, and `symdim` is the first dimension
begin Fourier transformed.

# Examples

## 1D, real to complex
```julia
A = JopFft(Float64,128)
m = rand(domain(A))
d = A*m
```
## 1D, complex to complex
```julia
A = JopFft(ComplexF64,128)
m = rand(domain(A))
d = A*m
```

## 2D, real to complex
```julia
A = JopFft(Float64,128,256)
m = rand(domain(A))
d = A*m
```

## 3D, real to complex, transforming over only the first dimension
```julia
A = JopFft(Float64,128,256; dims=(1,))
m = rand(domain(A))
d = A*m
```
"""
function JopFft(::Type{T},n::Vararg{Int,N};dims=()) where {T<:Real,N}
    dims = isempty(dims) ? ntuple(i->i, N) : dims
    JopLn(dom = JetSpace(T,n), rng = symspace(typeof(JopFft_df!),T,dims[1],n...), df! = JopFft_df!, df′! = JopFft_df′!, s = (dims=dims,))
end
function JopFft(::Type{T},n::Vararg{Int,N};dims=()) where {T<:Complex,N}
    dims = isempty(dims) ? ntuple(i->i, N) : dims
    JopLn(dom = JetSpace(T,n), rng = JetSpace(T,n), df! = JopFft_df!, df′! = JopFft_df′!, s = (dims=dims,))
end
JopFft(R::JetSpace; dims=()) = JopFft(eltype(R), size(R)...; dims=dims)
export JopFft

function JopFft_df!(d::Jets.SymmetricArray{Complex{T}}, m::AbstractArray{T}; kwargs...) where {T}
    JopFft_df!(parent(d), m; kwargs...)
    d
end
function JopFft_df!(d::Array{Complex{T}}, m::AbstractArray{T}; dims, kwargs...) where {T<:Real}
    P = plan_rfft(m, dims)
    JopFft_mul!(d, P, m, dims)
    d
end

function JopFft_df!(d::AbstractArray{T}, m::AbstractArray{T}; dims, kwargs...) where {T<:Complex}
    P = plan_fft(m, dims)
    JopFft_mul!(d, P, m, dims)
    d
end

function JopFft_mul!(d::AbstractArray{T}, P::FFTW.FFTWPlan, m::AbstractArray, dims) where {T}
    mul!(d, P, m)
    sc = fftscale(T, size(m), dims)
    d .*= sc
end

JopFft_df′!(m::AbstractArray{T}, d::Jets.SymmetricArray{Complex{T}}; kwargs...) where {T} = JopFft_df′!(m, parent(d); kwargs...)
function JopFft_df′!(m::AbstractArray{T}, d::Array{Complex{T}}; dims, kwargs...) where {T<:Real}
    P = plan_brfft(d, size(m, dims[1]), dims)
    JopFft_mul_adjoint!(m, P, copy(parent(d)), dims) # in-place version of brfft modifies the input buffer
    m
end

function JopFft_df′!(m::AbstractArray{T}, d::AbstractArray; dims, kwargs...) where {T<:Complex}
    P = plan_bfft(d, dims)
    JopFft_mul_adjoint!(m, P, d, dims)
    m
end

function JopFft_mul_adjoint!(m::AbstractArray{T}, P::FFTW.FFTWPlan, d::AbstractArray, dims) where {T}
    mul!(m, P, d)
    sc = fftscale(T, size(m), dims)
    m .*= sc
end

fftscale(::Type{T}, n, dims) where {T} = T(1.0/sqrt(mapreduce(dim->n[dim], *, dims)))

function _JopFft_symspace_map(idim, I, symdim, n)
    if idim == symdim
        return 2 + n[idim] - I[idim]
    else
        return I[idim]
    end
end

function JopFft_symspace_map(I, symdim::Int, n::NTuple{N,Int}, _n::NTuple{N,Int}) where {N}
    if I[symdim] > _n[symdim]
        J = ntuple(idim->_JopFft_symspace_map(idim, I, symdim, n), N)
        return CartesianIndex(J)
    else
        return CartesianIndex(I)
    end
end

function Jets.symspace(JopFft, ::Type{T}, symdim::Int, n::Vararg{Int,N}) where {T,N}
    _n = ntuple(idim->idim == symdim ? div(n[idim], 2) + 1 : n[idim], N)
    Jets.JetSSpace(Complex{T}, n, _n, I->JopFft_symspace_map(I, symdim, n, _n))
end


text:
###########################################################################
#  Semi-explicit Riemannian manifold Lagrangian Monte Carlo (RMLMC)
#
#  Reference: Lan S, Stathopoulos V, Shahbaba B, Girolami M. Lagrangian Dynamical Monte Carlo. arXiv, 2012
#
#  Parameters :
#    - nLeaps : number of leapfrog steps
#    - leapStep : leapfrog step size
#    - nNewton: number of Newton steps
#    - tuner: for tuning the RMLMC parameters
#
###########################################################################

export RMLMC

println("Loading RMLMC(nLeaps, leapStep, nNewton, tuner) sampler")

###########################################################################
# RMLMC specific tuners
###########################################################################
abstract RMLMCTuner <: MCMCTuner

###########################################################################
# RMLMC type
###########################################################################
immutable RMLMC <: MCMCSampler
  nLeaps::Int
  leapStep::Float64
  nNewton::Int
  tuner::Union(Nothing, MCMCTuner)

  function RMLMC(nLeaps::Int, leapStep::Float64, nNewton::Int, tuner::Union(Nothing, MCMCTuner))
    @assert nLeaps>0 "Number of leapfrog steps should be > 0"
    @assert leapStep>0 "Leapfrog step size should be > 0"
    @assert nNewton>0 "Number of Newton steps should be > 0"
    new(nLeaps, leapStep, nNewton, tuner)
  end
end

RMLMC(tuner::Union(Nothing, MCMCTuner)=nothing) = RMLMC(6, 0.5, 4, tuner)
RMLMC(nLeaps::Int, tuner::Union(Nothing, MCMCTuner)=nothing) = RMLMC(nLeaps, 3/nLeaps, 4, tuner)
RMLMC(nLeaps::Int, leapStep::Float64, tuner::Union(Nothing, MCMCTuner)=nothing) = RMLMC(nLeaps, leapStep, 4, tuner)
RMLMC(nLeaps::Int, nNewton::Int, tuner::Union(Nothing, MCMCTuner)=nothing) = RMLMC(nLeaps, 3/nLeaps, nNewton, tuner)
RMLMC(leapStep::Float64, tuner::Union(Nothing, MCMCTuner)=nothing) = RMLMC(int(floor(3/leapStep)), leapStep, 4, tuner)
RMLMC(leapStep::Float64, nNewton::Int, tuner::Union(Nothing, MCMCTuner)=nothing) =
  RMLMC(int(floor(3/leapStep)), leapStep, nNewton, tuner)

###########################################################################
# RMLMC sampler
###########################################################################
function SamplerTask(model::MCMCModel, sampler::RMLMC, runner::MCMCRunner)
  # Define local variables
  local pars, proposedPars
  local logTarget, proposedLogTarget
  local grad, proposedGrad
  local E, proposedE
  local G, proposedG, invG, proposedInvG, cholG, proposedCholG, dG, proposeddG
  local traceInvGxdG, proposedTraceInvGxdG
  local dphi, proposeddphi, C, proposedC
  local velocity, proposedVelocity, initLeapVelocity, vxC, deltaLogDet
  local nRandomLeaps
  local nLeaps, leapStep

  @assert hasgradient(model) "RMLMC sampler requires model with gradient function"
  @assert hastensor(model) "RMLMC sampler requires model with tensor function"
  @assert hasdtensor(model) "RMLMC sampler requires model with function of tensor derivatives"

   #  Task reset function
  function reset(resetPars::Vector{Float64})
    pars = copy(resetPars)
    logTarget, grad, G, dG = model.evalalldt(pars)
  end
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, reset)

  # Initialization
  pars = copy(model.init)
  logTarget, grad, G, dG = model.evalalldt(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

  invG = inv(G)
  cholG = chol(G)
  traceInvGxdG = Float64[trace(invG*dG[:, :, i]) for i = 1:model.size]
  dphi = -grad+0.5*traceInvGxdG
  C = 0.5*(permutedims(dG, [3 2 1])+permutedims(dG, [1 3 2])-dG)

  if isa(sampler.tuner, EmpMCTuner); tune = EmpiricalHMCTune(sampler.nLeaps, sampler.leapStep, 0, 0); end

  i = 1
  while true
    if isa(sampler.tuner, EmpMCTuner)
      tune.proposed += 1
      nLeaps, leapStep = tune.nLeaps, tune.leapStep
    else
      nLeaps, leapStep = sampler.nLeaps, sampler.leapStep
    end

    proposedPars = copy(pars)
    proposedG = copy(G)
    proposedInvG = copy(invG)
    proposedCholG = copy(cholG)
    proposeddG = copy(dG)
    proposedTraceInvGxdG = copy(traceInvGxdG)
    proposeddphi = copy(dphi)
    proposedC = copy(C)

    # Sample velocity and calculate current Energy
    velocity = (chol(proposedInvG))'*randn(model.size)
    proposedVelocity = copy(velocity)
    E = -logTarget+sum(log(diag(cholG)))+0.5*velocity'*G*velocity

    vxC = zeros(model.size, model.size)
    # Accumulate determinant to be adjusted in acceptance rate
    deltaLogDet = 0

    # Perform leapfrog steps
    nRandomLeaps = ceil(rand()*nLeaps)

    for j = 1:nRandomLeaps
      # Update velocity
      leapVelocity = copy(proposedVelocity)
      for k = 1:sampler.nNewton
        for r = 1:model.size
          vxC[r, :] = leapVelocity'*proposedC[:, :, r]
        end
        leapVelocity = proposedVelocity-(0.5*leapStep)*proposedInvG*(vxC*leapVelocity+proposeddphi)
      end
      proposedVelocity = copy(leapVelocity)

      #Update volume correction term
      deltaLogDet = deltaLogDet-logdet(proposedG+(leapStep)*vxC)

      # Update parameters
      proposedPars = proposedPars+leapStep*proposedVelocity

      proposedLogTarget, proposedGrad, proposedG, proposeddG = model.evalalldt(proposedPars)

      proposedInvG = inv(proposedG)
      proposedCholG = chol(proposedG)
      proposedTraceInvGxdG = Float64[trace(proposedInvG*proposeddG[:, :, k]) for k = 1:model.size]

      proposeddphi = -proposedGrad+0.5*proposedTraceInvGxdG
      proposedC = 0.5*(permutedims(proposeddG, [3 2 1])+permutedims(proposeddG, [1 3 2])-proposeddG)

      # Update velocity
      for k = 1:model.size
        vxC[k, :] = proposedVelocity'*proposedC[:, :, k]
      end
      deltaLogDet = deltaLogDet+logdet(proposedG-leapStep*vxC)

      proposedVelocity = proposedVelocity-(0.5*leapStep)*proposedInvG*(vxC*proposedVelocity+proposeddphi)
    end

    # Calculate energy based on the proposed parameters
    proposedE = -proposedLogTarget+sum(log(diag(proposedCholG)))+0.5*proposedVelocity'*proposedG*proposedVelocity

    # Accept according to ratio
    ratio = (E-proposedE)[1]+deltaLogDet

    if ratio > 0 || (ratio > log(rand()))  # i.e. if accepted
      produce(MCMCSample(proposedPars, proposedLogTarget, proposedGrad, pars, logTarget, grad, {"accept" => true}))
      pars, logTarget, grad = copy(proposedPars), copy(proposedLogTarget), copy(proposedGrad)
      G, dG, invG, cholG = copy(proposedG), copy(proposeddG), copy(proposedInvG), copy(proposedCholG)
      traceInvGxdG, C, dphi = copy(proposedTraceInvGxdG), copy(proposedC), copy(proposeddphi)
      if isa(sampler.tuner, EmpMCTuner); tune.accepted += 1; end
    else
      produce(MCMCSample(pars, logTarget, grad, pars, logTarget, grad, {"accept" => false}))
    end

    if isa(sampler.tuner, EmpMCTuner) && i<= runner.burnin && mod(i, sampler.tuner.adaptStep) == 0
      adapt!(tune, sampler.tuner)
      reset!(tune)
      if sampler.tuner.verbose
        println("Burn-in teration $i of $(runner.burnin): ", round(100*tune.rate, 2), " % acceptance rate")
      end
    end

    i += 1
  end
end


text:
@testset "Trivial tree" begin
    #= Manually create a simple tree with the following structure

        (1)
         |
        (2)
    =#
    leaf = BPLeaf("I'm a leaf", 1, :final)
    node = BPNode(0, [2])
    tree = BPTree(Dict(1 => node),
                  Dict(2 => leaf),
                  [2])

    # Check that printing does not error
    io = IOBuffer()
    println(io, leaf)
    println(io, node)
    println(io, tree)

    @test root(tree) == node
    @test nnodes(tree) == 2

    d = data(tree)
    @test length(d) == 1
    @test d[1] == "I'm a leaf"

    k = newid(tree)
    @test k ∉ 0:2  # The new ID must be new and not 0
    @test isa(k, Integer)

    discard_leaf!(tree, 2)
    # The last leaf was deleted so the tree should be empty
    @test nnodes(tree) == 0
end

@testset "Simple tree" begin
    #= Tests on a slighty less trivial tree

          (1)
         /   \
       (2)   (8)
      /  \     \
    (3)  (5)   (9)
     |   | \     \
    (4) (6) (7)  (10)
    =#

    n1 = BPNode(0, [2, 8])
    n2 = BPNode(1, [3, 5])
    n3 = BPNode(2, [4])
    n4 = BPLeaf("Leaf 4", 3, :working)
    n5 = BPNode(2, [6, 7])
    n6 = BPLeaf("Leaf 6", 5, :final)
    n7 = BPLeaf("Leaf 7", 5, :working)
    n8 = BPNode(1, [9])
    n9 = BPNode(8, [10])
    n10 = BPLeaf("Leaf 10", 9, :working)

    tree = BPTree(Dict(1 => n1, 2 => n2, 3 => n3, 5 => n5, 8 => n8, 9 => n9),
                  Dict(4 => n4, 6 => n6, 7 => n7, 10 => n10),
                  [3, 4, 8])

    # Check that printing does not error
    io = IOBuffer()
    println(io, tree)

    @test newid(tree) ∉ 0:10
    @test length(data(tree)) == 4

    @test nnodes(tree) == 10

    discard_leaf!(tree, 6)
    @test nnodes(tree) == 9

    discard_leaf!(tree, 4)
    @test nnodes(tree) == 7

    discard_leaf!(tree, 10)
    @test nnodes(tree) == 4

    @test length(data(tree)) == 1
end


text:
module TestColumnarArrayInterface

using FactCheck
using DataStructures

import Iterators: chain

import DBAPI
import DBAPI.ArrayInterfaces:
    ColumnarArrayInterface,
    ColumnarArrayQuery,
    ArrayInterfaceError

facts("Array interface") do
    # invalid
    @fact_throws ArrayInterfaceError Base.connect(ColumnarArrayInterface, [:foo], Vector[])
    @fact_throws ArrayInterfaceError Base.connect(identity, ColumnarArrayInterface, [:foo], Vector[])
    @fact_throws ArrayInterfaceError Base.connect(ColumnarArrayInterface, Symbol[], Vector[[1, 2, 3]])
    @fact_throws ArrayInterfaceError Base.connect(ColumnarArrayInterface, [:foo, :bar], Vector[[1], [2, 3]])

    # do block (without the do block)
    connection = Base.connect(identity, ColumnarArrayInterface, Symbol[], Vector[])
    @fact isa(connection, DBAPI.DatabaseConnection) --> true
    @fact DBAPI.isopen(connection) --> false

    # empty
    connection = Base.connect(ColumnarArrayInterface, Symbol[], Vector[])
    @fact isa(connection, DBAPI.DatabaseConnection) --> true
    @fact DBAPI.isopen(connection) --> true
    @fact DBAPI.commit(connection) --> exactly(nothing)
    @fact_throws DBAPI.NotSupportedError DBAPI.rollback(connection)
    @fact Base.isopen(connection) --> true

    cursor = DBAPI.cursor(connection)
    @fact isa(cursor, DBAPI.DatabaseCursor) --> true
    @fact isa(cursor, DBAPI.FixedLengthDatabaseCursor) --> true

    @fact_throws ArrayInterfaceError DBAPI.rows(cursor)
    @fact_throws ArrayInterfaceError DBAPI.columns(cursor)

    @fact_throws DBAPI.DatabaseQueryError DBAPI.execute!(cursor, ColumnarArrayQuery([:one], 1:0))
    @fact_throws DBAPI.DatabaseQueryError DBAPI.execute!(cursor, ColumnarArrayQuery(Symbol[], 1:1))
    @fact_throws DBAPI.DatabaseQueryError DBAPI.execute!(cursor, ColumnarArrayQuery([:one], 1:1))

    @fact DBAPI.execute!(cursor, ColumnarArrayQuery(Symbol[], 1:0)) --> exactly(nothing)

    @fact isempty(collect(DBAPI.rows(cursor))) --> true
    @fact isempty(collect(DBAPI.columns(cursor))) --> true
    @fact_throws BoundsError cursor[1, 1]
    @fact_throws BoundsError cursor[1, :one]
    @fact isempty(cursor) --> true
    @fact length(cursor) --> 0

    @fact_throws DBAPI.NotImplementedError cursor[:one, 1]
    @fact_throws DBAPI.NotImplementedError cursor[:one, :one]
    @fact_throws DBAPI.NotImplementedError cursor["one", "one"]

    @fact Base.close(connection) --> exactly(nothing)
    @fact Base.isopen(connection) --> false

    # non-empty
    connection = Base.connect(
        ColumnarArrayInterface,
        [:foo, :bar],
        Vector[[1, 2, 3], [3.0, 2.0, 1.0]]
    )
    @fact isa(connection, DBAPI.DatabaseConnection) --> true
    @fact DBAPI.isopen(connection) --> true
    @fact DBAPI.commit(connection) --> exactly(nothing)
    @fact_throws DBAPI.NotSupportedError DBAPI.rollback(connection)
    @fact Base.isopen(connection) --> true

    cursor = DBAPI.cursor(connection)
    @fact isa(cursor, DBAPI.DatabaseCursor) --> true
    @fact isa(cursor, DBAPI.FixedLengthDatabaseCursor) --> true

    @fact_throws DBAPI.DatabaseQueryError DBAPI.execute!(cursor, ColumnarArrayQuery([:one], 1:1))

    @fact DBAPI.execute!(cursor, ColumnarArrayQuery(Symbol[], 1:0)) --> exactly(nothing)

    @fact collect(DBAPI.rows(cursor)) --> isempty
    @fact collect(DBAPI.columns(cursor)) --> isempty
    @fact cursor --> isempty
    @fact length(cursor) --> 0
    @fact_throws BoundsError cursor[1, 1]
    @fact_throws BoundsError cursor[1, :one]

    @fact DBAPI.execute!(cursor, ColumnarArrayQuery([:foo, :bar], 1:3)) --> exactly(nothing)
    row_results = [
        (1, 3.0),
        (2, 2.0),
        (3, 1.0),
    ]
    nullable_row_results = map(x -> map(Nullable{Any}, x), row_results)
    @fact collect(DBAPI.rows(cursor)) --> row_results
    @fact collect(DBAPI.columns(cursor)) --> Vector[
        [1, 2, 3],
        [3.0, 2.0, 1.0],
    ]

    @fact cursor[3, :bar] --> 1.0
    @fact cursor[1, :foo] --> 1
    @fact_throws BoundsError cursor[4, :foo]
    @fact_throws BoundsError cursor[1, :one]
    @fact_throws BoundsError cursor[0, :one]

    @fact cursor --> not(isempty)
    @fact length(cursor) --> length(row_results)

    empty_data_structures = (
        Array{Any}[Array{Any}(0)],
        Array{Any}[Array{Any}(0, 0)],
        Dict{Any, Any}[Dict{Any, Any}()],
        PriorityQueue[PriorityQueue()],
        Dict{Any,Array{Any}}(1=>Array{Any}(0)),
        Dict{Any,Array{Any}}(1=>Array{Any}(0, 0)),
        Array{Nullable{Any}}[Array{Nullable{Any}}(0)],
        Array{Nullable{Any}}[Array{Nullable{Any}}(0, 0)],
        Dict{Any, Nullable{Any}}[Dict{Any, Nullable{Any}}()],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(0)),
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(0, 0)),
    )

    empty_2d_data_structures = (
        Array{Any}(0,0),
        Dict{Any, Any}(),
        PriorityQueue(),
        Array{Nullable{Any}}(0,0),
        Dict{Any, Nullable{Any}}(),
    )

    filled_pq = PriorityQueue()
    filled_pq[1] = 1

    filled_data_structures = (
        Array{Any}[Array{Any}(1)],
        Array{Any}[Array{Any}(1, 1)],
        Dict{Any, Any}[Dict{Any, Any}(1=>5)],
        PriorityQueue[filled_pq],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(1)),
        Array{Nullable{Any}}[Array{Nullable{Any}}(1)],
        Array{Nullable{Any}}[Array{Nullable{Any}}(1, 1)],
        Dict{Any, Nullable{Any}}[Dict{Any, Nullable{Any}}(1=>5)],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(1)),
    )

    filled_2d_pq = PriorityQueue()
    filled_2d_pq[1, 1] = 1

    filled_2d_data_structures = (
        Array{Any}(1),
        Array{Any}(1, 1),
        Dict{Any, Any}((1, 1)=>5),
        filled_2d_pq,
        Array{Nullable{Any}}(1),
        Array{Nullable{Any}}(1, 1),
        Dict{Any, Nullable{Any}}((1, 1)=>Nullable{Any}(5)),
    )

    for ds in empty_data_structures
        # ds will be a collection of an empty row
        @fact DBAPI.fetchintorows!(ds, cursor) --> (ds, 1)
    end

    for ds in empty_data_structures
        @fact DBAPI.fetchintocolumns!(ds, cursor) --> exactly((ds, 0))
    end

    for ds in empty_2d_data_structures
        @fact DBAPI.fetchintoarray!(ds, cursor) --> exactly((ds, 0))
    end

    for ds in filled_data_structures
        @fact DBAPI.fetchintorows!(ds, cursor) --> exactly((ds, 1))
    end

    for ds in filled_data_structures
        @fact DBAPI.fetchintocolumns!(ds, cursor) --> exactly((ds, 1))
    end

    for ds in filled_2d_data_structures
        @fact DBAPI.fetchintoarray!(ds, cursor) --> exactly((ds, 1))
    end

    first_empty(a::Associative) = isempty(first(values(a)))
    first_empty(a) = isempty(first(a))

    for empty_ds in empty_data_structures
        iters = 0
        for ds in DBAPI.DatabaseFetcher(:rows, empty_ds, cursor)
            # ds will be a collection containing an empty row
            @fact first_empty(ds) --> true
            iters += 1
        end
        @fact iters --> length(cursor)
    end

    for empty_ds in empty_data_structures
        for ds in DBAPI.DatabaseFetcher(:columns, empty_ds, cursor)
            @fact true --> false  # should never be reached
        end
    end

    for empty_ds in empty_2d_data_structures
        for ds in DBAPI.DatabaseFetcher(:array, empty_ds, cursor)
            @fact true --> false  # should never be reached
        end
    end

    comparison_data(::Any) = row_results
    comparison_data(::Nullable) = nullable_row_results

    for filled_ds in filled_data_structures
        for (idx, ds) in enumerate(DBAPI.DatabaseFetcher(:rows, filled_ds, cursor))
            item = ds[1][1]
            @fact item --> exactly(comparison_data(item)[idx][1])
        end
    end

    for filled_ds in filled_data_structures
        for (idx, ds) in enumerate(DBAPI.DatabaseFetcher(:columns, filled_ds, cursor))
            item = ds[1][1]
            @fact item --> exactly(comparison_data(item)[idx][1])
        end
    end

    for filled_ds in filled_2d_data_structures
        for (idx, ds) in enumerate(DBAPI.DatabaseFetcher(:array, filled_ds, cursor))
            item = ds[1, 1]
            @fact item --> exactly(comparison_data(item)[idx][1])
            @fact length(ds) --> 1
        end
    end

    # nullable db data
    connection = Base.connect(
        ColumnarArrayInterface,
        [:foo, :bar],
        Vector[[Nullable(1), Nullable(2), Nullable(3)], [3.0, 2.0, 1.0]]
    )
    @fact isa(connection, DBAPI.DatabaseConnection) --> true
    @fact DBAPI.isopen(connection) --> true
    @fact DBAPI.commit(connection) --> exactly(nothing)
    @fact_throws DBAPI.NotSupportedError DBAPI.rollback(connection)
    @fact Base.isopen(connection) --> true

    cursor = DBAPI.cursor(connection)
    @fact isa(cursor, DBAPI.DatabaseCursor) --> true
    @fact isa(cursor, DBAPI.FixedLengthDatabaseCursor) --> true

    @fact_throws DBAPI.DatabaseQueryError DBAPI.execute!(cursor, ColumnarArrayQuery([:one], 1:1))

    @fact DBAPI.execute!(cursor, ColumnarArrayQuery(Symbol[], 1:0)) --> exactly(nothing)

    @fact collect(DBAPI.rows(cursor)) --> isempty
    @fact collect(DBAPI.columns(cursor)) --> isempty
    @fact cursor --> isempty
    @fact length(cursor) --> 0
    @fact_throws BoundsError cursor[1, 1]
    @fact_throws BoundsError cursor[1, :one]

    @fact DBAPI.execute!(cursor, ColumnarArrayQuery([:foo, :bar], 2:3)) --> exactly(nothing)
    row_results = Any[
        (Nullable(2), 2.0),
        (Nullable(3), 1.0),
    ]
    nullable_row_results = map(x -> map(y -> isa(y, Nullable) ? y : Nullable{Any}(y), x), row_results)
    @fact isequal(collect(DBAPI.rows(cursor)), row_results) --> true
    @fact isequal(collect(DBAPI.columns(cursor)), Vector[
        [Nullable(2), Nullable(3)],
        [2.0, 1.0],
    ]) --> true

    # unsupported query types
    @fact_throws DBAPI.NotSupportedError DBAPI.execute!(
        cursor,
        "[:foo, :bar], 1:3",
    )
    @fact_throws DBAPI.NotSupportedError DBAPI.execute!(
        cursor,
        "[:foo, :bar], 1:3",
        (),
    )
    @fact_throws DBAPI.NotSupportedError DBAPI.execute!(
        cursor,
        DBAPI.StringMultiparameterQuery(
            "[:foo, :bar], 1:3",
            ((),())
        ),
    )
end

end


text:
using JLD, GibbsTypePriors, DataFramesMeta, RCall, Optim, OptimTestProblems, StatsBase, Optim, OptimTestProblems

include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/common_functions.jl")

R"library(tidyverse)
library(latex2exp)
library(viridis)
"

function approximation_prior_distribution(beta_,sigma,N,Nt)
    Pkn_numeric_ = GibbsTypePriors.Pkn_2PD.(1:N, N, beta_, sigma)
    Pkn_order2_ = GibbsTypePriors.Pkn_PY_pred_approx(N, beta_, sigma)
    Pkn_PY_Mult = GibbsTypePriors.Pkn_PYM.(1:N, N, Nt, beta_, sigma)
    R"SB_PY = Prior_on_K_SB_PY($beta_, $sigma, $N, $Nt, runs= 2*10^2)
    pk_sb_py = SB_PY$pk"
    @rget pk_sb_py
    df= DataFrame(Pkn_numeric = Pkn_numeric_,
                  Pkn_order2 = Pkn_order2_,
                  Pkn_PYM = Pkn_PY_Mult,
                  Pkn_SB  = pk_sb_py)
       return df
end

R"
StickBreakingPY <- function(alpha,sigma, N) {
  a<- rep(1-sigma,N-1)
  b<- alpha+ sigma*(1:N-1)
  V <- rbeta(N-1 , a, b)
  p    <- vector('numeric',length=N)
  p[1] <- V[1]
  for(l in 2:(N - 1))p[l] <- prod(1 - V[1:(l - 1)])*V[l]
  p[N] <- prod(1 - V)
  p
  return(p)
}

Prior_on_K_SB_PY<- function(alpha, sigma, N_s=100,N_tr=100, runs=10^4 ){
  array_nc<-c()
  i=1
  while (i<=runs){
    weights_PY<- StickBreakingPY(alpha,sigma,N_tr)
    c <- sample(1:N_tr,size=N_s, replace=TRUE, prob=weights_PY )
    n_c<- length(unique(c))
    array_nc[i]<- n_c
    i=i+1
  }
  p_k = tibble(k=as.numeric(names(table(array_nc))),
               p_k=as.vector(table(array_nc))/sum(table(array_nc)))
  p_zeros= tibble(k=(1:N_s)[!1:N_s%in%p_k$k],
                  p_k=rep(0,length((1:N_s)[!1:N_s%in%p_k$k])))
  pk_padded= rbind(p_k, p_zeros)%>% arrange(k)
  E_k=sum((1:N_s) *(pk_padded$p_k))
  V_k = sum((((1:N_s) - E_k)^2) *(pk_padded$p_k))
  return(list(pk=pk_padded$p_k, Ek= E_k,Vk=V_k ))
}
"

n=100
β= 1.0
ntr=250
sigma_vec= [0.25,0.75]
DF_all_100_py = map(x ->approximation_prior_distribution(β,x,n,ntr),sigma_vec)

n=100
β= 10.0
ntr=250
sigma_vec= [0.25,0.75]
DF_all_10_100_py = map(x ->approximation_prior_distribution(β,x,n,ntr),sigma_vec)

DF_100_250_py = [DF_all_100_py,DF_all_10_100_py]

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_5.jld", "DF_100_250_py" ,DF_100_250_py)

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
"""
    JuliaLikeInterpreter <: AbstractInterpreter

Abstract type for julia-like interpreter. The subtype of it usually
modifies the native julia interpreter a little bit by overloading
certain abstract interpretation interface, but forward most of the
interfaces to the native interpreter.
"""
abstract type JuliaLikeInterpreter <: AbstractInterpreter end

"""
    parent(interp::JuliaLikeInterpreter)

Return the native interpreter of Julia.
"""
Base.parent(interp::JuliaLikeInterpreter) = interp.native_interpreter

Core.Compiler.InferenceParams(interp::JuliaLikeInterpreter) = InferenceParams(interp.native_interpreter)
Core.Compiler.OptimizationParams(interp::JuliaLikeInterpreter) = OptimizationParams(interp.native_interpreter)
Core.Compiler.get_world_counter(interp::JuliaLikeInterpreter) = get_world_counter(interp.native_interpreter)
Core.Compiler.get_inference_cache(interp::JuliaLikeInterpreter) =
    get_inference_cache(interp.native_interpreter)
Core.Compiler.code_cache(interp::JuliaLikeInterpreter) = Core.Compiler.code_cache(interp.native_interpreter)
Core.Compiler.may_optimize(interp::JuliaLikeInterpreter) =
    Core.Compiler.may_optimize(interp.native_interpreter)
Core.Compiler.may_discard_trees(interp::JuliaLikeInterpreter) =
    Core.Compiler.may_discard_trees(interp.native_interpreter)
Core.Compiler.may_compress(interp::JuliaLikeInterpreter) =
    Core.Compiler.may_compress(interp.native_interpreter)
Core.Compiler.unlock_mi_inference(interp::JuliaLikeInterpreter, mi::Core.MethodInstance) =
    Core.Compiler.unlock_mi_inference(interp.native_interpreter, mi)
Core.Compiler.lock_mi_inference(interp::JuliaLikeInterpreter, mi::Core.MethodInstance) =
    Core.Compiler.lock_mi_inference(interp.native_interpreter, mi)
Core.Compiler.add_remark!(::JuliaLikeInterpreter, ::Core.Compiler.InferenceState, ::String) =
    nothing # println(msg)
@static if VERSION >= v"1.7.0-DEV.577"
    Core.Compiler.verbose_stmt_info(::JuliaLikeInterpreter) = false
end

function Core.Compiler.optimize(interp::JuliaLikeInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    @static if VERSION < v"1.8-"
        nargs = Int(opt.nargs) - 1
        ir = Core.Compiler.run_passes(opt.src, nargs, opt)
    else
        ir = Core.Compiler.run_passes(opt.src, opt)
    end

    ir = optimize(interp, opt, ir)
    @static if VERSION < v"1.7-DEV"
        Core.Compiler.finish(opt, params, ir, result)
    else
        Core.Compiler.finish(interp, opt, params, ir, result)
    end
end

"""
    optimize(interp::JuliaLikeInterpreter[, state::OptimizationState], ir::IRCode)

This method is for overloading, it will be executed after running Julia optimizers.
If you wish to customize the default Julia optimization passes, consider overloading
`Core.Compiler.optimize(interp, opt, params, result)`.
"""
function optimize end
optimize(interp::JuliaLikeInterpreter, ::OptimizationState, ir::IRCode) = optimize(interp, ir)
optimize(::JuliaLikeInterpreter, ir::IRCode) = ir


text:
"""
    tf!(dtm::AbstractMatrix{Real}, tf::AbstractMatrix{AbstractFloat})

Overwrite `tf` with the term frequency of the `dtm`.

Works correctly if `dtm` and `tf` are same matrix.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf!(dtm::AbstractMatrix{T1}, tf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    for i in 1:n
        words_in_document = 0
        for j in 1:p
            words_in_document += dtm[i, j]
        end
        tf[i, :] = dtm[i, :] ./ max(words_in_document, one(T1))
    end

    return tf
end

"""
    tf!(dtm::SparseMatrixCSC{Real}, tf::SparseMatrixCSC{AbstractFloat})

Overwrite `tf` with the term frequency of the `dtm`.

`tf` should have the has same nonzeros as `dtm`.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf!(dtm::SparseMatrixCSC{T}, tf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfvals = nonzeros(tf)
    @assert size(dtmvals) == size(tfvals)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = sum(dtm,dims=2)

    n, p = size(dtm)
    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfvals[j] = dtmvals[j] / max(words_in_documents[row], one(T))
       end
    end
    return tf
end

tf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf!(dtm, dtm)

tf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf!(dtm, dtm)

"""
    tf(dtm::DocumentTermMatrix)
    tf(dtm::SparseMatrixCSC{Real})
    tf(dtm::Matrix{Real})

Compute the `term-frequency` of the input.

# Example

```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
              StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)

julia> tf(m)
2×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 10 stored entries:
  [1, 1]  =  0.166667
  [2, 1]  =  0.166667
  [1, 2]  =  0.333333
  [2, 3]  =  0.333333
  [1, 4]  =  0.166667
  [2, 4]  =  0.166667
  [1, 5]  =  0.166667
  [2, 5]  =  0.166667
  [1, 6]  =  0.166667
  [2, 6]  =  0.166667
```

See also: [`tf!`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
tf(dtm::DocumentTermMatrix) = tf(dtm.dtm)

tf(dtm::Matrix{T}) where {T <: Real} = tf!(dtm, Array{Float64}(undef, size(dtm)...))

tf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf!(dtm, similar(dtm, Float64))

"""
    tf_idf!(dtm::AbstractMatrix{Real}, tf_idf::AbstractMatrix{AbstractFloat})

Overwrite `tf_idf` with the tf-idf (Term Frequency - Inverse Doc Frequency) of the `dtm`.

`dtm` and `tf-idf` must be matrices of same dimensions.

See also: [`tf`](@ref), [`tf!`](@ref) , [`tf_idf`](@ref)
"""
function tf_idf!(dtm::AbstractMatrix{T1}, tfidf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    tf!(dtm, tfidf)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    # TF-IDF is the product of TF and IDF
    for i in 1:n
        for j in 1:p
           tfidf[i, j] = tfidf[i, j] * idf[j]
        end
    end

    return tfidf
end

"""
    tf_idf!(dtm::SparseMatrixCSC{Real}, tfidf::SparseMatrixCSC{AbstractFloat})

Overwrite `tfidf` with the tf-idf (Term Frequency - Inverse Doc Frequency) of the `dtm`.

The arguments must have same number of nonzeros.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf_idf!(dtm::SparseMatrixCSC{T}, tfidf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfidfvals = nonzeros(tfidf)
    @assert size(dtmvals) == size(tfidfvals)

    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(dtm, dims=2))
    oneval = one(F)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfidfvals[j] = dtmvals[j] / max(words_in_documents[row], oneval) * idf[i]
       end
    end

    return tfidf
end

"""
    tf_idf!(dtm)

Compute tf-idf for `dtm`
"""
tf_idf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf_idf!(dtm, dtm)

tf_idf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf_idf!(dtm, dtm)

# This does not make sense, since DocumentTermMatrix is based on an array of integers
#tf_idf!(dtm::DocumentTermMatrix) = tf_idf!(dtm.dtm)

"""
    tf(dtm::DocumentTermMatrix)
    tf(dtm::SparseMatrixCSC{Real})
    tf(dtm::Matrix{Real})

Compute `tf-idf` value (Term Frequency - Inverse Document Frequency) for the input.

In many cases, raw word counts are not appropriate for use because:

- Some documents are longer than other documents
- Some words are more frequent than other words

A simple workaround this can be done by performing `TF-IDF` on a `DocumentTermMatrix`

# Example

```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
              StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)

julia> tf_idf(m)
2×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 10 stored entries:
  [1, 1]  =  0.0
  [2, 1]  =  0.0
  [1, 2]  =  0.231049
  [2, 3]  =  0.231049
  [1, 4]  =  0.0
  [2, 4]  =  0.0
  [1, 5]  =  0.0
  [2, 5]  =  0.0
  [1, 6]  =  0.0
  [2, 6]  =  0.0
```

See also: [`tf!`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
tf_idf(dtm::DocumentTermMatrix) = tf_idf(dtm.dtm)

tf_idf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf_idf!(dtm, similar(dtm, Float64))

tf_idf(dtm::Matrix{T}) where {T <: Real} = tf_idf!(dtm, Array{Float64}(undef, size(dtm)...))


text:
using Test
if @isdefined(EA_AS_PKG)
    import EscapeAnalysis: code_escapes, @code_escapes
    using EscapeAnalysis
else
    using Core.Compiler.EscapeAnalysis
    import Base: code_escapes
    import InteractiveUtils: @code_escapes
end
import Core: Argument, SSAValue

@static if isdefined(Core.Compiler, :alloc_array_ndims)
    import Core.Compiler: alloc_array_ndims
else
    function alloc_array_ndims(name::Symbol)
        if name === :jl_alloc_array_1d
            return 1
        elseif name === :jl_alloc_array_2d
            return 2
        elseif name === :jl_alloc_array_3d
            return 3
        elseif name === :jl_new_array
            return 0
        end
        return nothing
    end
end

isT(T) = (@nospecialize x) -> x === T
issubT(T) = (@nospecialize x) -> x <: T
isreturn(@nospecialize x) = isa(x, Core.ReturnNode) && isdefined(x, :val)
isthrow(@nospecialize x) = Meta.isexpr(x, :call) && Core.Compiler.is_throw_call(x)
isnew(@nospecialize x) = Meta.isexpr(x, :new)
isϕ(@nospecialize x) = isa(x, Core.PhiNode)
function with_normalized_name(@nospecialize(f), @nospecialize(x))
    if Meta.isexpr(x, :foreigncall)
        name = x.args[1]
        nn = EscapeAnalysis.normalize(name)
        return isa(nn, Symbol) && f(nn)
    end
    return false
end
isarrayalloc(@nospecialize x) = with_normalized_name(nn->!isnothing(alloc_array_ndims(nn)), x)
isarrayresize(@nospecialize x) = with_normalized_name(nn->!isnothing(EscapeAnalysis.array_resize_info(nn)), x)
isarraycopy(@nospecialize x) = with_normalized_name(nn->EscapeAnalysis.is_array_copy(nn), x)
import Core.Compiler: argextype, singleton_type
const EMPTY_SPTYPES = Any[]
iscall(y) = @nospecialize(x) -> iscall(y, x)
function iscall((ir, f), @nospecialize(x))
    return iscall(x) do @nospecialize x
        Core.Compiler.singleton_type(Core.Compiler.argextype(x, ir, EMPTY_SPTYPES)) === f
    end
end
iscall(pred::Function, @nospecialize(x)) = Meta.isexpr(x, :call) && pred(x.args[1])

let setup_ex = quote
        mutable struct SafeRef{T}
            x::T
        end
        Base.getindex(s::SafeRef) = getfield(s, 1)
        Base.setindex!(s::SafeRef, x) = setfield!(s, 1, x)

        mutable struct SafeRefs{S,T}
            x1::S
            x2::T
        end
        Base.getindex(s::SafeRefs, idx::Int) = getfield(s, idx)
        Base.setindex!(s::SafeRefs, x, idx::Int) = setfield!(s, idx, x)
    end
    global function EATModule(setup_ex = setup_ex)
        M = Module()
        Core.eval(M, setup_ex)
        return M
    end
    Core.eval(@__MODULE__, setup_ex)
end


text:
function unit_value(parameter::T) where T <: NumberParameter
    (parameter.value - parameter.min) / (parameter.max - parameter.min)
end

function unit_value(parameter::T) where T <: EnumParameter
    (parameter.value - 1) / (length(parameter.values) - 1)
end

function unit_value!(parameter::T, unit::Float64) where T <: EnumParameter
    @assert 0.0 <= unit <= 1.0
    value = round((unit * (length(parameter.values) - 1)) + 1)
    @inbounds parameter.current = parameter.values[parameter.value]
    parameter.value
end

function unit_value!(parameter::IntegerParameter, unit::Float64)
    @assert 0.0 <= unit <= 1.0
    value = round(Int, (unit * (parameter.max - parameter.min)) + parameter.min)
    parameter.value
end

function unit_value!(parameter::FloatParameter, unit::Float64)
    @assert 0.0 <= unit <= 1.0
    value = (unit * (parameter.max - parameter.min)) + parameter.min
    parameter.value
end


text:
using Lumberjack

export Wfst, Arc, add_arc!, add_initial_state!, add_final_state!

# The arc type that specify "arcs", "edges" or "transition rules".
type Arc
    from
    to
    input::String
    output::String
    weight::Float64
end

# The weighted FST type.
type Wfst
    states::Set
    input_alphabet::Set{String}
    output_alphabet::Set{String}
    initial_states::Set
    final_states::Set
    # Arcs of the form: (from, to, input, output, weight)
    arcs::Set{Arc}
    # Maps from initial and final states to their weights.
    initial_weights::Dict
    final_weights::Dict
    semiring::Semiring
end

# Empty constructor
Wfst(semiring::Semiring) = Wfst(Set(), Set{String}(), Set{String}(), Set(),
        Set(), Set{Arc}(), Dict(), Dict(), semiring)

# Add a state to the initial states list with the given weight.
function add_initial_state!(wfst::Wfst, state, weight::Float64)
    @assert state in wfst.states
    wfst.initial_states = union(wfst.initial_states, Set([state]))
    wfst.initial_weights[state] = weight
end

# Add a state to the final states list with the given weight.
function add_final_state!(wfst::Wfst, state, weight::Float64)
    @assert state in wfst.states
    wfst.final_states = union(wfst.final_states, Set([state]))
    wfst.final_weights[state] = weight
end

function add_arc!(wfst::Wfst,
        from, to,
        input::String, output::String, weight::Float64)
    wfst.states = union(wfst.states, Set([[from], [to]]))
    if !ismatch(r"<.*>", input)
        wfst.input_alphabet = union(wfst.input_alphabet, Set([input]))
    end
    if !ismatch(r"<.*>", output)
        wfst.output_alphabet = union(wfst.output_alphabet, Set([output]))
    end
    wfst.arcs =
            union(wfst.arcs, Set([Arc(from, to, input, output, weight)]))
end


text:
using Printf
using Plots

using Oceananigans
using Oceananigans.Advection
using JULES

using JULES: IdealGas
using Oceananigans.Grids: Cell, xnodes

ENV["GKSwstype"] = "100"

N = 512
L = 1
T = 1e-1
Δt = 2e-4
Nt = Int(T/Δt)

ideal_gas = IdealGas(Float64, JULES.R⁰, 1; T₀=1, p₀=1, s₀=1, u₀=0)

topo = (Periodic, Periodic, Periodic)
domain = (x=(0, L), y=(0, 1), z=(0, 1))
grid = RegularCartesianGrid(topology=topo, size=(N, 1, 1), halo=(4, 4, 4); domain...)

model = CompressibleModel(
                      grid = grid,
                 advection = WENO5(),
                     gases = (ρ=ideal_gas,),
    thermodynamic_variable = Energy(),
                   closure = IsotropicDiffusivity(ν=0, κ=0),
                   gravity = 0.0
)

g  = model.gravity
gas = model.gases.ρ
R, cₚ, cᵥ = gas.R, gas.cₚ, gas.cᵥ
u₀₀, T₀₀, ρ₀₀, s₀₀ = gas.u₀, gas.T₀, gas.ρ₀, gas.s₀

ρ₀(x, y, z) = 1
p₀(x, y, z) = 1
T₀(x, y, z) = p₀(x, y, z) / (R*ρ₀(x, y, z))
ρe₀(x, y, z) = ρ₀(x, y, z) * (u₀₀ + cᵥ * (T₀(x, y, z) - T₀₀) + g*z)

@inline x′(x, t, L, U) = mod(x + L/2 - U * t, L) - L/2
@inline ϕ_Gaussian(x, t; L, U, a=1, c=1/8) = a * exp(-x′(x, t, L, U)^2 / (2c^2))
@inline ϕ_Square(x, t; L, U, w=0.15) = -w <= x′(x, t, L, U) <= w ? 1.0 : 0.0
@inline ϕ_Sine(x, t; L, U) = sin(2π * x)
ρu₀(x, y, z) = ϕ_Sine(x, 0, L=L, U=1)

set!(model.tracers.ρ, ρ₀)
set!(model.tracers.ρe, ρe₀)
set!(model.momenta.ρu, ρu₀)
update_total_density!(model)

anim = @animate for n in 1:Nt
    @info "iteration $n/$Nt"
    time_step!(model, Δt)

    title = @sprintf("Burgers equation t=%.3f", model.clock.time)

    x = xnodes(Cell, grid)
    ρu = interior(model.momenta.ρu)[:]
    plot(x, ρu, lw=2, label="", title=title, xlims=(0, L), ylims=(-2, 2), dpi=200)
end

mp4(anim, "burgers_equation.mp4", fps=60)


text:
## Code for our Types and conversion methods

## Symbol class for controlling dispatch
abstract SymbolicObject <: Number

## Basic types defined here
"""

* `Sym("x")`
* `Sym(:x)`
* `Sym("x, y, z")`

The `Sym` type is an immutable type holding a reference to an
underlying python-based SymPy object. Many methods are extended to the
`Sym` type. Instances can be constructed in many ways. The one caveat
is the variables can not be function names in base.

"""
immutable Sym <: SymbolicObject
    x::PyCall.PyObject
end
Sym(s::SymbolicObject) = s

## Matrix type
immutable SymMatrix <: SymbolicObject
    x::PyCall.PyObject
end

## sets
immutable SymSet <: SymbolicObject
    x::PyCall.PyObject
end

## complex float
## this cause issue with printing on non-complex objects
#mpctype = sympy.mpmath["ctx_mp_python"]
#pytype_mapping(mpctype, Sym)

## some typealiases
typealias SymOrReal @compat Union{Sym, Real}
typealias SymOrNumber @compat Union{Sym, Number}
typealias SymOrString @compat Union{Sym, AbstractString}
typealias SymbolicTypes  @compat Union{AbstractString, Symbol, SymbolicObject}

PyCall.PyObject(x::SymbolicObject) = x.x

## Promotion
## promote up to symbolic so that math ops work
promote_rule{T<:SymbolicObject, S<:Number}(::Type{T}, ::Type{S} ) = T

## Conversion
convert(::Type{Sym}, o::PyCall.PyObject) = Sym(o)
convert(::Type{PyObject}, s::Sym) = s.x

function convert(::Type{Tuple}, o::PyCall.PyObject)
    ## check that o is a tuple?
    ## PyCall.pytypeof(o)
    n = o[:__len__]()
    ntuple(i -> o[:__getitem__](i-1), n)
end

## rational
convert{T<:SymbolicObject}(::Type{T}, x::Rational) = sympy_meth(:Rational, x.num, x.den)::T

## big. Need mpmath installed separately -- not as a SymPy module as that is how it is called in PyCall
convert{T<:SymbolicObject}(::Type{T}, x::BigFloat) = Sym(PyCall.PyObject(x))::T
convert(::Type{Sym}, x::Complex{BigFloat}) = Sym(PyCall.PyObject(x))::Sym

## real
convert{S<:SymbolicObject, T <: Real}(::Type{S}, x::T) = sympy_meth(:sympify, x)::S
convert{T <: Real}(::Type{T}, x::Sym) = convert(T, project(x))

## complex
## IM is SymPy's "i" (sympy[:I], not Python's
## Sym(PyCall.PyObject(im)) which gives 1j.
function convert(::Type{Sym}, x::Complex)
    y = ifelse(isa(x, Complex{Bool}), real(x) + imag(x) * im, x)
    real(y) + imag(y) * IM
end
convert(::Type{Complex}, x::Sym) = complex(map(x -> convert(Float64, x), x[:as_real_imag]())...)::Sym
complex(x::Sym) = convert(Complex, x)
complex(xs::Array{Sym}) = map(complex, xs)

## matrices
convert(::Type{SymMatrix}, o::PyCall.PyObject) = SymMatrix(o)
convert(::Type{Sym}, o::SymMatrix) = Sym(o.x)
convert(::Type{SymMatrix}, o::Sym) = SymMatrix(o.x)

## string
convert(::Type{Sym}, o::AbstractString) = sympy_meth(:sympify, o)
convert(::Type{Sym}, o::Symbol) = sympy_meth(:sympify, string(o))

"""

Get the free symbols in a more convenient form than as returned by `free_symbols`.

Just `free_symbols` now does the same thing. This function will be deprecated.
"""
function get_free_symbols(ex::Sym)
    free = free_symbols(ex)
    n = free[:__len__]()
    n == 0 && return(Sym[])

    vars = [free[:pop]()]
    if n > 1
        for i in 1:(n-1)
            push!(vars, free[:pop]())
        end
    end
    vars
end

## Conversion to function a bit hacky
## we use free_symbols to get the free symbols, then create a function
## with arguments in this order. No problem with only one variable, but
## may be confusing when more than one in ex.
## Output is symbolic. Conversion is necessary to use output as Julia values.
## SymPy has `lamdify` for this task too.
function convert(::Type{Function}, ex::Sym)
    vars = free_symbols(ex)
    len = length(vars)
    if len == 0
        return x -> ex
    end
    local out
    (args...) -> begin
        out = ex
        for i in 1:len
            out = object_meth(out, :subs, vars[i], args[i]) #convert(Function, out[:subs])(project(vars[i]), args[i])
        end
        out
    end
end

## For plotting we need to know if a function has 1, 2, or 3 free variables
function as_nfunction(ex::Sym, nvars=1)
    free = free_symbols(ex)
    vars = [free[:pop]()]
    for i in 1:free[:__len__]()
        push!(vars, free[:pop]())
    end
    len = length(vars)

    if len == nvars
        convert(Function, ex)
    else
        throw(DimensionMismatch("Expecting $nvars free variables and found $len"))
    end
end

"""
Convert to a julia-valued function from R -> R::Float64
"""
type ScalarFunction; end
function Base.convert(::Type{ScalarFunction}, ex::Sym)
    vars = free_symbols(ex)
    length(vars) == 1 || error("Scalar function conversion is for expression of one variable")
    u -> convert(Float64, N(subs(ex, vars[1], u)))
end


text:
# Must be run in a separate process from the rest of the tests!

using RuntimeGeneratedFunctions
using Serialization

RuntimeGeneratedFunctions.init(@__MODULE__)

f = @RuntimeGeneratedFunction(:(x->"Hi from a separate process. x=$x"))

serialize(stdout, f)


text:
using Distributions
using CausalInference
using LightGraphs
using DelimitedFiles
using LinearAlgebra

if !@isdefined(pathNCI)
    pathNCI = download("http://nugget.unisa.edu.au/ParallelPC/data/real/NCI-60.csv")
end
data = readdlm(pathNCI, ',')
d = size(data,2)
X = data[:,1:d]'
d, n = size(X)
C = Symmetric(cor(X, dims=2))
p = 0.01
@time h, s = skeleton(d, gausscitest, (C,n), quantile(Normal(), 1-p/2))
println("inferred edges ", ne(h))


text:
"""
    badge(id, name, department)

Return the employee badge label.
"""
function badge(id, name, department)
    department = isnothing(department) ? "owner" : department
    prefix = isnothing(id) ? "" : "[$id] - "

    prefix * "$name - $(uppercase(department))"
end


text:
module ParTemp

import Base.append!
using Distributions

function metropolisMV(curr, logprob, S)
  cand = rand(MvNormal(curr, S))
  logU = log(rand())
  logP = logprob(cand) - logprob(curr)
  accept = logP > logU
  draw = accept ? cand : curr
  return (draw, accept)
end

mutable struct Buffer{T}
  n::Int
  x::Vector{T}
  Buffer(T::Type, n::Int) = new{T}(n, T[])
end

tail(x, n) = n > length(x) ? x : x[(end-n+1):end]

function append!(buffer::Buffer, x)
  append!(buffer.x, x)
  buffer.x = tail(buffer.x, buffer.n)
end

end


text:
struct SphereSignedDistanceFunction{T<:Real} <: AbstractSignedDistanceFunction
    radius::T
end

radius(sdf::SphereSignedDistanceFunction) = sdf.radius
(sdf::SphereSignedDistanceFunction{T1})(point::Point{3, T2}) where {T1<:Real, T2<:Real} =
    norm(point) - T2(radius(sdf))

# Currently always pointing out
normal(::SphereSignedDistanceFunction{T1}, point::Point{3, T2}) where {T1<:Real, T2<:Real} =
    Vec3{T2}(normalize(point))