text:
import MechanicalSketch: @import_expand, color_with_lumin, empty_figure, finish
import MechanicalSketch: rect, squircle, PALETTE, ∙, WI, HE, O, EM
import MechanicalSketch: @latexify, latexify, place_image, modify_latex
import MechanicalSketch: paint, setmesh, mesh, box, rule, BoundingBox
import MechanicalSketch: @layer, color_with_lumin, setline, sethue, Length
import MechanicalSketch: scale_to_pt, pixelwidth, text, draw_background
import MechanicalSketch: Point, row_height, circle, blend, addstop, setblend
import MechanicalSketch: setdash, line, RGBA, chroma, saturation, color_with_saturation
import MechanicalSketch: settext, brush, translate, pixelwidth, box
import MechanicalSketch: setopacity, scale_to_pt, scale_pt_to_unit
import Base.pop!

if !@isdefined m²
    @import_expand(~m, s, °)
end
include("test_functions_36.jl")

empty_figure(filename = joinpath(@__DIR__, "test_36.png"), height = 15m);
pt = O + (-WI / 2, -HE / 2) + (EM, 2.5EM)
str = "<big>Testing</big>    a) latex    b) brush strokes    c) sequence of cropped panels\r
       The blocks use the golden ratio:"
settext(str, pt, markup = true)
pt += (0, 0.5EM)
# This assigns value to w:
expression = :(w = (1.0m * (1 + √5)) / 2)
eval(expression)
showexpression = :($expression = $(round(m, w; digits =3)))
ptul, ptbr, latexscale = place_image(pt + (3EM, 0), latexify(showexpression); height = 3EM, centered = false)
pt += (0, 1.5EM + ptbr[2] - ptul[2])
settext("Each 'comic panel', too. Panel origo is circled:", pt, markup = true)
# A panel has origo at the middle of the bottom line (that's how it will be drawn)
pt += (0.0m, - 2.5m)
panelorigo = pt + (3EM, 0) + (w, 0.0m)
circle(panelorigo; d= 1m)
pa = Container{Panel}(name = "panel")
moveto!(pa, panelorigo )
drawit(pa)
pt += (0, 2EM)
settext("Start comic strip:", pt, markup = true)
#
pa.name = "1"
move_y!(pa, -7m)
set_width_height!(pa, 3.5w, 3.5m)
moveto_x!(pa, 1.2 * 3.5w / 2 + scale_pt_to_unit(m) * (-WI / 2 ))
Δy = height(pa) * (2 /  (1 + √5))^3
ground = Container{Ground}(; y = Δy, w = 5 *  height(pa), h = Δy)
putin!(pa, ground)
stack = Container{Stack}(; x = -1m, name = "Stack")
putin!(ground, stack)
putin!(stack, Thing{Section}{Mounted}(;name = "1st"))
drawit(pa)

pa.name = "2"
move_x!(pa,  1.2 * 3.5w)
putin!(ground, Thing{Section}{Transit}(;name = "2nd", y = 0.25m, x = 1.0m ))
drawit(pa)

#
pa.name = "3"
move_x!(pa,  1.2 * 3.5w)
putin!(stack, Thing{Section}{Mounted}(;name = "2nd"))
pop!(ground)
drawit(pa)

finish()

#=import MechanicalSketch: BezierPath, BezierPathSegment, bezier, bezier′, bezier′′, makebezierpath, drawbezierpath, bezierpathtopoly, beziertopoly, pathtobezierpaths,
bezierfrompoints, beziercurvature, bezierstroke, setbezierhandles, shiftbezierhandles
h = 1m
empty_figure(filename = joinpath(@__DIR__, "test_36.png"), height = 15m);
pts = [Point(-w,-w),
Point(-w,w),
Point(w,w),
Point(w,-w),
Point(-w,-w)]

for (sm, ptx) in zip(range(0.0, 1.0; length = 4), 0.5.*[-0.75WI, -0.25WI, 0.25WI, 0.75WI])
   pt = Point(ptx, 0.0)
   be = makebezierpath(pts .+ pt; smoothing = sm)
   drawbezierpath(be, :stroke; close = true)
end
be = makebezierpath(pts; smoothing = 0.1)
=#


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

pde = Array{Any, 1}(undef, 2);
mesh = Array{Any, 1}(undef, 2);

# create pde and mesh for each PDE model
include("pdeapp1.jl");
include("pdeapp2.jl");

# call exasim to generate and run C++ code to solve the PDE models
sol,pde,mesh,master,dmd,compilerstr,runstr = exasim(pde,mesh);

# visualize the numerical solution of the PDE model using Paraview
for m = 1:length(pde)
    pde[m].visscalars = ["temperature", 1];  # list of scalar fields for visualization
    pde[m].visvectors = ["temperature gradient", [2, 3]]; # list of vector fields for visualization
    pde[m].visfilename = "dataout" * string(m) * "/output";
    vis(sol[m],pde[m],mesh[m]); # visualize the numerical solution
end
print("Done!");

ext: jl


text:
struct FINITO_lbfgs_adaptive_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg, TH, Ts} <: CIAO_iterable
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    η::R                    # factor for updating gammas
    sweeping::Int8             # to only use one stepsize γ
    batch::Int                # batch size
    α::R                    # in (0, 1), e.g.: 0.99
    H::TH
    adaptive::Bool          # adaptive stepsizes
    tol_b::R
    S::Ts            # functions s_i for calculating the Hassian of f_i s
end

mutable struct FINITO_lbfgs_adaptive_state{R<:Real,Tx, TH}
    γ::Tx                    # stepsize parameter
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    H::TH                   # Hessian approx
    # av_x::Tx                # sum x_i
    sum_fx::R
    sum_fv::R
    sum_x::Tx               # sum x_i/γ_i
    sum_nrmx::R               # sum |x_i|^2/γ_i
    sum_nabla::Tx            # sum \nabla f(x_i)
    sum_nabla_fz::Tx            # sum \nabla f(x_i)
    sum_innprod::R           # sum  < \nabla f(x_i), x_i >
    f_x::R                  # sum f(x_i)
    # some extra placeholders
    val_fg::R                   # value of g
    z::Tx
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
    zbar::Tx                # bar z
    zbar_prev::Maybe{Tx}    # bar z previous
    vbar::Tx            # v
    vbar_prev::Maybe{Tx}# v
    dir::Tx                 # direction
    # sum_nabla::Tx              # for linesearch
    z_trial::Tx             # linesearch candidate
    inds::Array{Int}        # needed for shuffled only!
    τ::Float64              # number of epochs
end

function FINITO_lbfgs_adaptive_state(γ::Tx, hat_γ::R, av::Tx, ind, d, H::TH, sum_x, sum_nrmx, sum_nabla, sum_innprod, f_x) where {R,Tx,TH}
    return FINITO_lbfgs_adaptive_state{R,Tx,TH}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        H,
        R(0),
        R(0),
        sum_x,
        sum_nrmx,
        sum_nabla,
        copy(av),
        sum_innprod,
        f_x,
        R(0),
        copy(av),
        copy(av),
        copy(av),
        copy(av),
        nothing,
        copy(av),
        nothing,
        copy(av),
        copy(av),
        collect(1:d),
        1.0
        )
end

function Base.iterate(iter::FINITO_lbfgs_adaptive_iterable{R}) where {R}
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))

    #initializing the vectors
    sum_nabla = zero(iter.x0)
    sum_innprod = R(0)
    f_x = R(0)
    for i = 1:N  # nabla f(x0)
        ∇f, fi_x = gradient(iter.F[i], iter.x0)
        sum_innprod += real(dot(∇f, iter.x0))
        sum_nabla .+= ∇f
        f_x += fi_x
    end

    # updating the stepsize
    if iter.γ === nothing ##  Lip. of f(x) = 1/N sum nabla f_i(x)
        if iter.L === nothing ||  iter.adaptive
           xeps = iter.x0 .+ one(R)
           av_eps = zero(iter.x0)
           for i in 1:N
                ∇f, ~ = gradient(iter.F[i], xeps)
                av_eps .+= ∇f
            end
            nmg = norm(sum_nabla - av_eps)
            t = 1
            while nmg < eps(R)  # in case xeps has the same gradient
                println("initial upper bound for L is too small")
                xeps = iter.x0 .+ rand(t * [-1, 1], size(iter.x0))
                av_eps = zero(iter.x0)
                for i in 1:N
                    ∇f, ~ = gradient(iter.F[i], xeps)
                    av_eps .+= ∇f
                end
                # grad_f_xeps, f_xeps = gradient(iter.F[i], xeps)
                nmg = norm(sum_nabla - av_eps)
                t *= 2
            end
            ###### decide how to initialize!
            L_int = nmg / (t * sqrt(length(iter.x0)))
            L_int /= iter.N # to account for 1/N
            γ = iter.α / (L_int)
            # γ = iter.S(iter.x0)
            println("gamma specified by L_int: ", γ)
        else
            println("gamma specified by L")
            isa(iter.L, R) ? (γ = iter.α * R(iter.N) / iter.L) :
                (γ = iter.α * R(N) / (maximum(iter.L)) )
        end
    else
        @warn "stepsize γ unspecified"
    end
    γ = fill(γ, (N,))
    hat_γ = 1 / sum(1 ./ γ)
    println("gggggg: ", hat_γ)
    # update the average
    av = copy(sum_nabla) .* (-hat_γ / N)
    av .+= iter.x0

    sum_x = iter.x0 /hat_γ
    sum_nrmx = N * norm(iter.x0)^2
    state = FINITO_lbfgs_adaptive_state(γ, hat_γ, av, ind, cld(N, r), iter.H, sum_x, sum_nrmx, sum_nabla, sum_innprod, f_x)

    return state, state
end

function Base.iterate(
    iter::FINITO_lbfgs_adaptive_iterable{R},
    state::FINITO_lbfgs_adaptive_state{R},
) where {R}

    if state.zbar_prev === nothing
        state.zbar_prev = zero(state.z)
        state.vbar_prev = zero(state.z)
    end

    # ------------- reset gamma by hessian initialization -------------
    # γ_prev = state.γ
    # state.γ = iter.S(state.z)
    # state.av .-=  ((state.γ - γ_prev) / iter.N) .* state.sum_nabla

    # full z update
     while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 1, ($(state.γ))"
            return nothing
        end
        prox!(state.zbar, iter.g, state.av, state.hat_γ)

        f_model_cnst = iter.N * state.sum_nrmx / 2 + state.f_x - state.sum_innprod

        f_model_z = iter.N * real(dot(state.zbar, state.sum_x))
        f_model_z -=  iter.N * norm(state.zbar)^2 / (2 * state.hat_γ)
        f_model_z -= real(dot(state.zbar, state.sum_nabla))
        for i in 1:iter.N
            f_model_z += iter.F[i](state.zbar)
        end

        tol = 10^(-6)  * (1 + abs(f_model_cnst))
        # tol = eps(R)
        R(f_model_z) <= f_model_cnst + tol && break
        println("ls 1")
        println("x: ", f_model_z)
        println("y: ", f_model_cnst)
        println("γ: ", state.γ)
        println("sum_nrmx: ", state.sum_nrmx)

        γ_prev = state.hat_γ
        state.γ *= iter.η
        state.hat_γ *= iter.η
        # updates
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla
        state.sum_nrmx ./= iter.η
        state.sum_x ./= iter.η
        reset!(state.H)
    end

    state.f_x = R(0)
    state.sum_innprod = R(0)
    state.sum_nabla = zero(state.av)
    for i = 1:iter.N
        state.∇f_temp, fi_z = gradient(iter.F[i], state.zbar) # update the gradient
        state.f_x += fi_z
        state.sum_nabla .+= state.∇f_temp
    end
    state.sum_innprod += real(dot(state.sum_nabla,state.zbar))
    state.av .= state.zbar .- (state.hat_γ / iter.N) .* state.sum_nabla
    state.sum_nrmx = norm(state.zbar)^2 * iter.N
    state.sum_nabla_fz .= state.sum_nabla

    # full v update
    while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 2, ($(state.γ))"
            return nothing
        end
        state.val_fg = prox!(state.vbar, iter.g, state.av, state.hat_γ) # bar v update
        state.sum_fv = 0
        for i in 1:iter.N
            state.sum_fv += iter.F[i](state.vbar)
        end

        f_model_cnst = iter.N * norm(state.vbar - state.zbar)^2/(2 * state.hat_γ)

        f_model_z = - real(dot(state.sum_nabla, state.vbar - state.zbar))
        f_model_z += state.sum_fv
        f_model_z -=  state.f_x

        tol = 10^(-6)  * (1 + abs(f_model_cnst))
        R(f_model_z) <= f_model_cnst + tol && break
        println("ls 2")

        γ_prev = state.hat_γ
        state.hat_γ *= 0.8
        state.γ *= 0.8
        # update av
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla_fz
        reset!(state.H)
    end
    state.vbar .-= state.zbar # \bar v- \bar z

    # prepare for linesearch: compute varphi(bar z)
    envVal = state.f_x / iter.N
    envVal += state.val_fg
    envVal += real(dot(state.sum_nabla, state.vbar)) / iter.N
    envVal += norm(state.vbar)^2 / (2 *  state.hat_γ)

    # update lbfgs
    update!(state.H, state.zbar - state.zbar_prev, -state.vbar +  state.vbar_prev)
    mul!(state.dir, state.H, state.vbar)
    # store vectors for next update
    copyto!(state.zbar_prev, state.zbar)
    copyto!(state.vbar_prev, state.vbar)

    state.τ = 1.0
    hat_γ_prev = state.hat_γ
    γ_prev = state.γ
    # while true  ### z_line should be removed using some other placeholder
    for i=1:5
        state.hat_γ = hat_γ_prev
        state.γ = γ_prev
        state.z_trial .=  state.zbar .+ (1- state.τ) .* state.vbar + state.τ * state.dir

        # compute varphi(z_trial)
        state.sum_nabla = zero(state.temp)
        state.sum_fx = 0
        for i = 1:iter.N
            state.∇f_temp, fi_z = gradient(iter.F[i], state.z_trial) # update the gradient
            state.sum_fx += fi_z
            state.sum_nabla .+= state.∇f_temp
        end

        while true
            if state.hat_γ < iter.tol_b / iter.N
                @warn "parameter `γ` became too small at ls 3, ($(state.γ))"
                return nothing
            end
            state.av .= state.z_trial
            state.av .-= (state.hat_γ / iter.N) .* state.sum_nabla
            state.val_fg = prox!(state.z, iter.g, state.av, state.hat_γ)

            temp = 0
            for i = 1:iter.N
                ~, fi_z = gradient(iter.F[i], state.z) # update the gradient
                temp += fi_z
            end

            envVal_trial = 0
            envVal_trial += state.sum_fx
            state.z .-= state.z_trial # r_γ(z^trial)
            envVal_trial += real(dot(state.sum_nabla, state.z))
            envVal_trial += iter.N * norm(state.z)^2 / (2 *  state.hat_γ)

            tol = 10^(-6)  * (1 + abs(envVal_trial))
            if temp <= envVal_trial + tol
                break
            end
            state.hat_γ *= 0.8
            state.γ *= 0.8
            println("ls 3")
        end
        envVal_trial /= iter.N
        envVal_trial += state.val_fg

        envVal_trial <= envVal + eps(R) && break
        state.τ /= 50       ##### bug prone: change in reporting if you change this! ######
        println("ls on τ")
    end
    state.zbar .= state.z_trial # of line search
    ###### at the momoent I am computing the next prox twice (it is already computed in the ls)
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled

    state.f_x = R(0)
    state.sum_innprod =  R(0)
    state.sum_nrmx =  R(0)
    state.sum_x = zero(state.temp)

    for j in state.inds
        for i in state.ind[j]
            # backtrack γ (warn if γ gets too small)

            hat_γ_prev = state.hat_γ
            γ_prev = state.γ[i]
            state.γ[i] = iter.S(state.zbar, i)

            while true
                state.hat_γ = 1/sum(1/state.γ)

                if state.hat_γ < iter.tol_b / iter.N
                    @warn "parameter `γ` became too small ($(state.γ))"
                    return nothing
                end

                println("norm av: ", norm(state.av))
                state.av .= state.av + (hat_γ_prev - state.hat_γ)/iter.N .* state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.zbar
                prox!(state.z, iter.g, state.av, state.hat_γ)

                fi_x = gradient!(state.∇f_temp, iter.F[i], state.zbar) # update the gradient
                state.val_fg = iter.F[i](state.z)  ##### change to function value! #####
                state.temp .= state.z .- state.zbar # placeholder for z - x^full

                fi_model = fi_x + real(dot(state.∇f_temp, state.temp)) +
                    (0.5 * iter.α * iter.N / state.γ[i]) * (norm(state.temp)^2)

                tol = 10^(-6)  * (1 + abs(fi_model))
                R(state.val_fg) <= fi_model + tol && break

                println("ls 4: ", state.hat_γ)
                println("x: ", state.val_fg)
                println("y: ", fi_model)
                println("sum_nabla: ", norm(state.sum_nabla))
                println("norm: ", state.γ[i])
                println("norm: ", norm(state.av))
                println("norm: ", norm(state.zbar))
                println("av+: ", norm((hat_γ_prev - state.hat_γ)/iter.N .* state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.zbar))
                println("dot: ", (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev))
                println("dot: ", (hat_γ_prev - state.hat_γ)/iter.N)
                return nothing
                γ_prev = state.γ[i]
                hat_γ_prev = state.hat_γ

                state.γ[i] *= iter.η
                reset!(state.H)
            end

            # iterates
            state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
            state.sum_nabla .-= state.∇f_temp # update sum nabla f_i for next iter
            gradient!(state.∇f_temp, iter.F[i], state.z) # update the gradient
            state.sum_nabla .+= state.∇f_temp # update sum nabla f_i for next iter
            state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
            state.av .+= state.hat_γ * (state.z .- state.zbar) ./ state.γ[i]

            # updates for the next linesearch
            state.f_x += state.val_fg
            state.sum_x += state.z / state.γ[i]
            state.sum_nrmx   += norm(state.z)^2 /state.γ[i]
            state.sum_innprod += real(dot(state.∇f_temp, state.z))
        end
    end
    # println("current stepsize is $(state.γ)")
    return state, state
end
solution(state::FINITO_lbfgs_adaptive_state) = state.z
epoch_count(state::FINITO_lbfgs_adaptive_state) = state.τ   # number of epochs is 2+ 1/tau , where 1/tau is from ls

#### TODO:
    # careful with the minibatch and adaptive stepsizes!
    # stepsizes if gamma is given...
    # add LS parameter as iter field


text:
struct MaternAcv{D,L} <: TimeSeriesModel{D,Float64}
    σ::SHermitianCompact{D,Float64,L}
    ν::SHermitianCompact{D,Float64,L}
    a::SHermitianCompact{D,Float64,L}
    a²::SHermitianCompact{D,Float64,L}
    νplushalf::SHermitianCompact{D,Float64,L}
    variance_part::SHermitianCompact{D,Float64,L}
    sdfconst::SHermitianCompact{D,Float64,L}
    acvconst::SHermitianCompact{D,Float64,L}
    function MaternAcv{D,L}(θ) where {D,L}
        @boundscheck checkparameterlength(θ, MaternAcv{D,L})
        L == triangularnumber(D) || error("MaternAcv{D,L} should satisfy L == D*(D+1)÷2")
        all(x->x>0, θ) || throw(ArgumentError("all parameters of MaternAcv should be > 0."))
        σ = @views SHermitianCompact(SVector{L,Float64}(θ[1:L]))
        all(i==j ? true : σ[i,j] < 1 for i in 1:size(σ,1) for j in 1:i) || throw(ArgumentError("ρ parameters must be < 1."))
        σ = @views SHermitianCompact(SVector{L,Float64}(θ[1:L]))
        ν = @views SHermitianCompact(SVector{L,Float64}(θ[L+1:2L]))
        a = @views SHermitianCompact(SVector{L,Float64}(θ[2L+1:end]))
        a² = SHermitianCompact((a.lowertriangle).^2)
        νplushalf = SHermitianCompact((ν.lowertriangle).+0.5)
        variance_part = SHermitianCompact(SMatrix{D,D,Float64}(i==j ? σ[i,i]^2 : σ[i,i]*σ[j,j]*σ[i,j] for i in 1:D, j in 1:D))
        sdfconst = SHermitianCompact(matern_sdf_normalising.(ν.lowertriangle,a.lowertriangle).*variance_part.lowertriangle)
        acvconst = SHermitianCompact(matern_acv_normalising.(ν.lowertriangle).*variance_part.lowertriangle)
        new{D,L}(σ,ν,a,a²,νplushalf,variance_part,sdfconst,acvconst)
    end
end
matern_acv_normalising(ν) = 2^(1-ν) / gamma(ν)

npars(::Type{MaternAcv{D,L}}) where {D,L} = 3triangularnumber(D)
lowerbounds(::Type{MaternAcv{D,L}}) where {D,L} = lowerbounds(Matern{D,L})
upperbounds(::Type{MaternAcv{D,L}}) where {D,L} = upperbounds(Matern{D,L})

parameternames(::Type{MaternAcv{D,L}}) where {D,L} = parameternames(Matern{D,L})

function add_sdf!(out, model::MaternAcv{D,L}, ω) where {D,L}
    ω² = ω^2
    count = 1
    for i ∈ 1:D, j ∈ i:D
        out[count] += model.sdfconst[i,j]/((model.a²[i,j]+ω²)^(model.νplushalf[i,j]))
        count += 1
    end
    return nothing
end

function acv!(out, model::MaternAcv{D,L}, τ::Number) where {D,L}
    modτ = abs(τ)
    count = 1
    for i ∈ 1:D, j ∈ i:D
        amodτ = model.a[i,j]*modτ
        out[count] = abs(τ)>1e-10 ? model.acvconst[i,j] * (amodτ)^model.ν[i,j] * besselk(model.ν[i,j],amodτ) : model.variance_part[i,j]
        count += 1
    end
    return nothing
end

struct MaternAcv1D <: TimeSeriesModel{1,Float64}
    σ::Float64
    ν::Float64
    a::Float64
    σ²::Float64
    a²::Float64
    νplushalf::Float64
    sdfconst::Float64
    acvconst::Float64
    function MaternAcv1D(σ,ν,a)
        σ > 0 || throw(ArgumentError("MaternAcv1D process requires 0 < σ."))
        ν > 0 || throw(ArgumentError("MaternAcv1D process requires 0 < ν."))
        a > 0 || throw(ArgumentError("MaternAcv1D process requires 0 < a."))
        sdfconst = matern_sdf_normalising(ν,a)*σ^2
        acvconst = matern_acv_normalising(ν)*σ^2
        new(σ, ν, a, σ^2, a^2, ν+0.5, sdfconst, acvconst)
    end
    function MaternAcv1D(x::AbstractVector{Float64})
        length(x) == npars(MaternAcv1D) || throw(ArgumentError("MaternAcv1D process has $(npars(Matern1D)) parameters, but $(length(x)) were provided."))
        @inbounds MaternAcv1D(x[1], x[2], x[3])
    end
end

npars(::Type{MaternAcv1D}) = 3

sdf(model::MaternAcv1D, ω) = model.sdfconst/((model.a²+ω^2)^(model.νplushalf))

acv(model::MaternAcv1D, τ) = abs(τ)>1e-10 ? model.acvconst * (amodτ)^model.ν * besselk(model.ν,a*abs(τ)) : model.σ²


text:
#miscellaneous helper functions and constants
Q95 = sqrt(2.)*erfinv(0.95)

range(x) = (minimum(x), maximum(x))

#95% normal quantile

roundv(v::Vector, i) = map(i -> round(i, 2), v)
roundv(v , i) = round(v,i)

eps2 = sqrt(eps())

intervalgaussian(U, a, b) = broadcast(intervalgaussian, U, a, b)
function intervalgaussian (U::Real,a::Real, b::Real)
	a, b = min(a,b), max(a,b)
	c =  erf(a/sqrt(2.))
	d =  erf(b/sqrt(2.))
	sqrt(2.)*erfinv(c + U.*(d-c))
end

function scalar(x)
assert(length(x) == 1)
x[1]
end

function norma!(x)
 minx = min(x)
 x = (x - minx) / (max(x)-minx)
 x
end

function cut(x, a, b)
	warn("cut(): better use clamp()")
	clamp(x, a,b)
end

#compute normal confidence interval for Monte Carlo estimate, precision aware rounding

function mc(M)
  m = mean(M)
  ste = std(M)/sqrt(length(M))
  round([m, Q95*ste], int(2-log(10,Q95*ste)))
end

#for sequential estimates, it suffices to provide the running sum Z, running sum of squares Z2 and the number of observations k

va(k, Z, Z2) =  Z2/k-(Z/k)^2
ste(k, Z, Z2) = sqrt(Z2/k^2-(Z/k)^2/k)

function mc2(k, Z, Z2, normalci = true)
	m = Z/k

	if normalci F = Q95 else F = 1. end

	try
		ste = sqrt(Z2/k-m.^2)/sqrt(k)

		res = (m, F*ste)

		return res
	catch
		return (m, NaN.*m)
	end

end

function selfadjmc(k, X, X2)
	m = X[1]/X[2]
	v = (va(k, X[1], X2[1,1]) - 2m*(X2[1,2]/k - X[1]*X[2]/k^2) + m^2*va(k, X[2], X2[2,2]))/k
	try
		ste = sqrt(v)
		res = [m, ste]
		return res
	catch
		return [m, NaN.*m]
	end

end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const OBJC_NEW_PROPERTIES = 1

# Skipping MacroDefinition: USER_ADDR_NULL ( ( user_addr_t ) 0 )
# Skipping MacroDefinition: CAST_USER_ADDR_T ( a_ptr ) ( ( user_addr_t ) ( ( uintptr_t ) ( a_ptr ) ) )

# const LITTLE_ENDIAN = __DARWIN_LITTLE_ENDIAN
# const BIG_ENDIAN = __DARWIN_BIG_ENDIAN
# const PDP_ENDIAN = __DARWIN_PDP_ENDIAN
# const BYTE_ORDER = __DARWIN_BYTE_ORDER

# Skipping MacroDefinition: ntohs ( x ) __DARWIN_OSSwapInt16 ( x )
# Skipping MacroDefinition: htons ( x ) __DARWIN_OSSwapInt16 ( x )
# Skipping MacroDefinition: ntohl ( x ) __DARWIN_OSSwapInt32 ( x )
# Skipping MacroDefinition: htonl ( x ) __DARWIN_OSSwapInt32 ( x )
# Skipping MacroDefinition: ntohll ( x ) __DARWIN_OSSwapInt64 ( x )
# Skipping MacroDefinition: htonll ( x ) __DARWIN_OSSwapInt64 ( x )
# Skipping MacroDefinition: NTOHL ( x ) ( x ) = ntohl ( ( __uint32_t ) x )
# Skipping MacroDefinition: NTOHS ( x ) ( x ) = ntohs ( ( __uint16_t ) x )
# Skipping MacroDefinition: NTOHLL ( x ) ( x ) = ntohll ( ( __uint64_t ) x )
# Skipping MacroDefinition: HTONL ( x ) ( x ) = htonl ( ( __uint32_t ) x )
# Skipping MacroDefinition: HTONS ( x ) ( x ) = htons ( ( __uint16_t ) x )
# Skipping MacroDefinition: HTONLL ( x ) ( x ) = htonll ( ( __uint64_t ) x )
# Skipping MacroDefinition: major ( x ) ( ( int32_t ) ( ( ( u_int32_t ) ( x ) >> 24 ) & 0xff ) )
# Skipping MacroDefinition: minor ( x ) ( ( int32_t ) ( ( x ) & 0xffffff ) )
# Skipping MacroDefinition: makedev ( x , y ) ( ( dev_t ) ( ( ( x ) << 24 ) | ( y ) ) )

# const NBBY = __DARWIN_NBBY
# const NFDBITS = __DARWIN_NFDBITS

# Skipping MacroDefinition: howmany ( x , y ) __DARWIN_howmany ( x , y )

# const FD_SETSIZE = __DARWIN_FD_SETSIZE

# Skipping MacroDefinition: FD_SET ( n , p ) __DARWIN_FD_SET ( n , p )
# Skipping MacroDefinition: FD_CLR ( n , p ) __DARWIN_FD_CLR ( n , p )
# Skipping MacroDefinition: FD_ZERO ( p ) __DARWIN_FD_ZERO ( p )
# Skipping MacroDefinition: FD_ISSET ( n , p ) __DARWIN_FD_ISSET ( n , p )
# Skipping MacroDefinition: FD_COPY ( f , t ) __DARWIN_FD_COPY ( f , t )

# const NULL = __DARWIN_NULL
# const RENAME_SECLUDE = 0x00000001
# const RENAME_SWAP = 0x00000002
# const RENAME_EXCL = 0x00000004
# const BUFSIZ = 1024
# const EOF = -1
# const FOPEN_MAX = 20
# const FILENAME_MAX = 1024
# const P_tmpdir = "/var/tmp/"
# const L_tmpnam = 1024
# const TMP_MAX = 308915776
# const SEEK_SET = 0
# const SEEK_CUR = 1
# const SEEK_END = 2
# const stdin = __stdinp
# const stdout = __stdoutp
# const stderr = __stderrp
# const L_ctermid = 1024

# Skipping MacroDefinition: getc_unlocked ( fp ) __sgetc ( fp )
# Skipping MacroDefinition: putc_unlocked ( x , fp ) __sputc ( x , fp )
# Skipping MacroDefinition: getchar_unlocked ( ) getc_unlocked ( stdin )
# Skipping MacroDefinition: putchar_unlocked ( x ) putc_unlocked ( x , stdout )
# Skipping MacroDefinition: fropen ( cookie , fn ) funopen ( cookie , fn , 0 , 0 , 0 )
# Skipping MacroDefinition: fwopen ( cookie , fn ) funopen ( cookie , 0 , fn , 0 , 0 )
# Skipping MacroDefinition: feof_unlocked ( p ) __sfeof ( p )
# Skipping MacroDefinition: ferror_unlocked ( p ) __sferror ( p )
# Skipping MacroDefinition: clearerr_unlocked ( p ) __sclearerr ( p )
# Skipping MacroDefinition: fileno_unlocked ( p ) __sfileno ( p )
# Skipping MacroDefinition: sprintf ( str , ... ) __builtin___sprintf_chk ( str , 0 , __darwin_obsz ( str ) , __VA_ARGS__ )
# Skipping MacroDefinition: snprintf ( str , len , ... ) __builtin___snprintf_chk ( str , len , 0 , __darwin_obsz ( str ) , __VA_ARGS__ )
# Skipping MacroDefinition: vsprintf ( str , format , ap ) __builtin___vsprintf_chk ( str , 0 , __darwin_obsz ( str ) , format , ap )
# Skipping MacroDefinition: vsnprintf ( str , len , format , ap ) __builtin___vsnprintf_chk ( str , len , 0 , __darwin_obsz ( str ) , format , ap )

const INT8_MAX = 127
const INT16_MAX = 32767
const INT32_MAX = 2147483647
const INT64_MAX = Int64(9223372036854775807)
const INT8_MIN = -128
const INT16_MIN = -32768
const INT32_MIN = -INT32_MAX - 1
const INT64_MIN = -INT64_MAX - 1
const UINT8_MAX = 255
const UINT16_MAX = 65535
const UINT32_MAX = UInt32(4294967295)
const UINT64_MAX = UInt64(@int128_str("18446744073709551615"))
const INT_LEAST8_MIN = INT8_MIN
const INT_LEAST16_MIN = INT16_MIN
const INT_LEAST32_MIN = INT32_MIN
const INT_LEAST64_MIN = INT64_MIN
const INT_LEAST8_MAX = INT8_MAX
const INT_LEAST16_MAX = INT16_MAX
const INT_LEAST32_MAX = INT32_MAX
const INT_LEAST64_MAX = INT64_MAX
const UINT_LEAST8_MAX = UINT8_MAX
const UINT_LEAST16_MAX = UINT16_MAX
const UINT_LEAST32_MAX = UINT32_MAX
const UINT_LEAST64_MAX = UINT64_MAX
const INT_FAST8_MIN = INT8_MIN
const INT_FAST16_MIN = INT16_MIN
const INT_FAST32_MIN = INT32_MIN
const INT_FAST64_MIN = INT64_MIN
const INT_FAST8_MAX = INT8_MAX
const INT_FAST16_MAX = INT16_MAX
const INT_FAST32_MAX = INT32_MAX
const INT_FAST64_MAX = INT64_MAX
const UINT_FAST8_MAX = UINT8_MAX
const UINT_FAST16_MAX = UINT16_MAX
const UINT_FAST32_MAX = UINT32_MAX
const UINT_FAST64_MAX = UINT64_MAX
const INTPTR_MIN = INT64_MIN
const INTPTR_MAX = INT64_MAX
const UINTPTR_MAX = UINT64_MAX
const INTMAX_MIN = INT64_MIN
const INTMAX_MAX = INT64_MAX
const UINTMAX_MAX = UINT64_MAX
const PTRDIFF_MIN = INT64_MIN
const PTRDIFF_MAX = INT64_MAX
const SIZE_MAX = UINT64_MAX
const RSIZE_MAX = SIZE_MAX >> 1
# const WCHAR_MAX = __WCHAR_MAX__
# const WCHAR_MIN = -WCHAR_MAX - 1
const WINT_MIN = INT32_MIN
const WINT_MAX = INT32_MAX
const SIG_ATOMIC_MIN = INT32_MIN
const SIG_ATOMIC_MAX = INT32_MAX

# Skipping MacroDefinition: INT8_C ( v ) ( v )
# Skipping MacroDefinition: INT16_C ( v ) ( v )
# Skipping MacroDefinition: INT32_C ( v ) ( v )
# Skipping MacroDefinition: INT64_C ( v ) ( v ## LL )
# Skipping MacroDefinition: UINT8_C ( v ) ( v ## U )
# Skipping MacroDefinition: UINT16_C ( v ) ( v ## U )
# Skipping MacroDefinition: UINT32_C ( v ) ( v ## U )
# Skipping MacroDefinition: UINT64_C ( v ) ( v ## ULL )
# Skipping MacroDefinition: INTMAX_C ( v ) ( v ## L )
# Skipping MacroDefinition: UINTMAX_C ( v ) ( v ## UL )

const CLOCKS_PER_SEC = 1000000

# begin enum ANONYMOUS_1
const ANONYMOUS_1 = UInt32
const _CLOCK_REALTIME = (UInt32)(0)
const _CLOCK_MONOTONIC = (UInt32)(6)
const _CLOCK_MONOTONIC_RAW = (UInt32)(4)
const _CLOCK_MONOTONIC_RAW_APPROX = (UInt32)(5)
const _CLOCK_UPTIME_RAW = (UInt32)(8)
const _CLOCK_UPTIME_RAW_APPROX = (UInt32)(9)
const _CLOCK_PROCESS_CPUTIME_ID = (UInt32)(12)
const _CLOCK_THREAD_CPUTIME_ID = (UInt32)(16)
# end enum ANONYMOUS_1

const CLOCK_REALTIME = _CLOCK_REALTIME
const CLOCK_MONOTONIC = _CLOCK_MONOTONIC
const CLOCK_MONOTONIC_RAW = _CLOCK_MONOTONIC_RAW
const CLOCK_MONOTONIC_RAW_APPROX = _CLOCK_MONOTONIC_RAW_APPROX
const CLOCK_UPTIME_RAW = _CLOCK_UPTIME_RAW
const CLOCK_UPTIME_RAW_APPROX = _CLOCK_UPTIME_RAW_APPROX
const CLOCK_PROCESS_CPUTIME_ID = _CLOCK_PROCESS_CPUTIME_ID
const CLOCK_THREAD_CPUTIME_ID = _CLOCK_THREAD_CPUTIME_ID
const RS_MAX_STRONG_SUM_LENGTH = 32
const RS_DEFAULT_BLOCK_LEN = 2048
const int8_t = UInt8
const int16_t = Int16
const int32_t = Cint
const int64_t = Clonglong
const u_int8_t = Cuchar
const u_int16_t = UInt16
const u_int32_t = UInt32
const u_int64_t = Culonglong
const register_t = Int64
# const intptr_t = __darwin_intptr_t
const uintptr_t = Culong
const user_addr_t = u_int64_t
const user_size_t = u_int64_t
const user_ssize_t = Int64
const user_long_t = Int64
const user_ulong_t = u_int64_t
const user_time_t = Int64
# NOTE: definition below is based only on above.
# TODO: audit Julia codebase to find actual definition
const time_t = user_time_t
const user_off_t = Int64
const syscall_arg_t = u_int64_t

const uint8_t = Cuchar
const uint16_t = UInt16
const uint32_t = UInt32
const uint64_t = Culonglong
const int_least8_t = Int8
const int_least16_t = Int16
const int_least32_t = Int32
const int_least64_t = Int64
const uint_least8_t = UInt8
const uint_least16_t = UInt16
const uint_least32_t = UInt32
const uint_least64_t = UInt64
const int_fast8_t = Int8
const int_fast16_t = Int16
const int_fast32_t = Int32
const int_fast64_t = Int64
const uint_fast8_t = UInt8
const uint_fast16_t = UInt16
const uint_fast32_t = UInt32
const uint_fast64_t = UInt64
const intmax_t = Clong
const uintmax_t = Culong

# mutable struct timespec
#     tv_sec::__darwin_time_t
#     tv_nsec::Clong
# end

mutable struct tm
    tm_sec::Cint
    tm_min::Cint
    tm_hour::Cint
    tm_mday::Cint
    tm_mon::Cint
    tm_year::Cint
    tm_wday::Cint
    tm_yday::Cint
    tm_isdst::Cint
    tm_gmtoff::Clong
    tm_zone::Cstring
end

const clockid_t = Void
const rs_byte_t = UInt8
const rs_long_t = intmax_t

# begin enum ANONYMOUS_2
const ANONYMOUS_2 = UInt32
const RS_DELTA_MAGIC = (UInt32)(1920139830)
const RS_MD4_SIG_MAGIC = (UInt32)(1920139574)
const RS_BLAKE2_SIG_MAGIC = (UInt32)(1920139575)
# end enum ANONYMOUS_2

const rs_magic_number = Void

# begin enum ANONYMOUS_3
const ANONYMOUS_3 = UInt32
const RS_LOG_EMERG = (UInt32)(0)
const RS_LOG_ALERT = (UInt32)(1)
const RS_LOG_CRIT = (UInt32)(2)
const RS_LOG_ERR = (UInt32)(3)
const RS_LOG_WARNING = (UInt32)(4)
const RS_LOG_NOTICE = (UInt32)(5)
const RS_LOG_INFO = (UInt32)(6)
const RS_LOG_DEBUG = (UInt32)(7)
# end enum ANONYMOUS_3

const rs_loglevel = Void

# begin enum rs_result
const rs_result = UInt32
const RS_DONE = (UInt32)(0)
const RS_BLOCKED = (UInt32)(1)
const RS_RUNNING = (UInt32)(2)
const RS_TEST_SKIPPED = (UInt32)(77)
const RS_IO_ERROR = (UInt32)(100)
const RS_SYNTAX_ERROR = (UInt32)(101)
const RS_MEM_ERROR = (UInt32)(102)
const RS_INPUT_ENDED = (UInt32)(103)
const RS_BAD_MAGIC = (UInt32)(104)
const RS_UNIMPLEMENTED = (UInt32)(105)
const RS_CORRUPT = (UInt32)(106)
const RS_INTERNAL_ERROR = (UInt32)(107)
const RS_PARAM_ERROR = (UInt32)(108)
# end enum rs_result

struct rs_stats
    op::Cstring
    lit_cmds::Cint
    lit_bytes::rs_long_t
    lit_cmdbytes::rs_long_t
    copy_cmds::rs_long_t
    copy_bytes::rs_long_t
    copy_cmdbytes::rs_long_t
    sig_cmds::rs_long_t
    sig_bytes::rs_long_t
    false_matches::Cint
    sig_blocks::rs_long_t
    block_len::Csize_t
    in_bytes::rs_long_t
    out_bytes::rs_long_t
    start::time_t
    _end::time_t
end

const rs_stats_t = Void

mutable struct rs_mdfour
end

const rs_mdfour_t = Void
const rs_weak_sum_t = UInt32
const rs_strong_sum_t = NTuple{32, Cuchar}

mutable struct rs_signature
end

const rs_signature_t = Void

mutable struct rs_buffers_s
    next_in::Cstring
    avail_in::Csize_t
    eof_in::Cint
    next_out::Cstring
    avail_out::Csize_t
end

const rs_buffers_t = Void

mutable struct rs_job
end

const rs_job_t = Void


text:
@test begin
    c = [Segment(Point(0, 0, 0), Point(1, 1, 1))]

    counter = BoxCounter(c)
    sizes, counts = [], []

    state = start(counter, 8)

    rate = Inf
    while !done(counter, state)
        state = next(counter, state)
        size, count = last(state)
        push!(sizes, size)
        push!(counts, count)

        if length(sizes) > 1
            rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
        end
        println((size, count, rate))
    end
    # abs(rate - 1) < 1E-2
    true
end

@test begin
    a = collect(linspace(0, 2, 10000))
    c = [Segment(Point(a[i], a[i], a[i]), Point(a[i+1], a[i+1], a[i+1]))
         for i in 1:length(a)-1]

    counter = BoxCounter(c)
    sizes, counts = [], []

    state = start(counter, 8)

    rate = Inf
    while !done(counter, state)
        state = next(counter, state)
        size, count = last(state)
        push!(sizes, size)
        push!(counts, count)

        if length(sizes) > 1
            rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
        end
        println((size, count, rate))
    end
    # abs(rate - 1) < 1E-2
    true
end

# # Check here that bbox counting gets approx the right dims
include("curves.jl")

@test begin
    θ = linspace(0, 2*pi, 10000)
    x = sin.(θ)
    y = cos.(θ)
    z = 0.5 + 0.5*x
    # Lift
    c = [Segment(Point(x[i], y[i], z[i]), Point(x[i+1], y[i+1], z[i+1])) for i in 1:10000-1]

    counter = BoxCounter(c)
    sizes, counts = [], []

    state = start(counter, 13)

    rate = Inf
    while !done(counter, state)
        state = next(counter, state)
        size, count = last(state)
        push!(sizes, size)
        push!(counts, count)

        if length(sizes) > 1
            rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
        end
        println((size, count, rate))
    end
    # abs(rate - 1) < 1E-2
    true
end

# # Diagonal of one piece
# @test begin
#     c = [Segment(Point(0, 0), Point(1, 1))]

#     counter = BoxCounter(c)
#     sizes, counts = [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#         end
#         # println((size, count, rate))
#     end
#     abs(rate - 1) < 1E-3
# end

# # Diagonal of more pieces
# @test begin
#     x = collect(linspace(0, 1, 10000))
#     c = [Segment(Point(x[i], x[i]), Point(x[i+1], x[i+1])) for i in 1:length(x)-1]

#     counter = BoxCounter(c)
#     sizes, counts = [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#         end
#         # println((size, count, rate))
#     end
#     abs(rate - 1) < 1E-3
# end

# # Circle
# @test begin
#     c = unit_circle(10000)

#     counter = BoxCounter(c)
#     sizes, counts = [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#         end
#         # println((size, count, rate))
#     end
#     abs(rate - 1) < 1E-3
# end

# # Koch
# @test begin
#     c = koch_flake(9)

#     counter = BoxCounter(c)
#     sizes, counts, rates = [], [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#             push!(rates, rate)
#         end
#         # println((size, count, rate))
#     end
#     approx = mean(rates[end-5:end])
#     exact = log(4)/log(3)
#     @show (approx, exact)
#     abs(approx - exact) < 1E-2
# end

# @test begin
#     c = koch_quadratic_1(8)

#     counter = BoxCounter(c)
#     sizes, counts, rates = [], [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#             push!(rates, rate)
#         end
#         # println((size, count, rate))
#     end
#     approx = mean(rates[end-5:end])
#     exact = log(5)/log(3)
#     @show (approx, exact)
#     abs(approx - exact) < 1E-2
# end

# @test begin
#     c = koch_quadratic_2(7)

#     counter = BoxCounter(c)
#     sizes, counts, rates = [], [], []

#     state = start(counter, 10)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#             push!(rates, rate)
#         end
#         # println((size, count, rate))
#     end
#     approx = round(mean(rates[end-3:end]), 2)
#     exact = 1.5
#     @show (approx, exact)
#     approx == exact
# end

# @test begin
#     c = dragon_curve(17)

#     counter = BoxCounter(c)
#     sizes, counts, rates = [], [], []

#     state = start(counter, 13)

#     rate = Inf
#     while !done(counter, state)
#         state = next(counter, state)
#         size, count = last(state)
#         push!(sizes, size)
#         push!(counts, count)

#         if length(sizes) > 1
#             rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
#             push!(rates, rate)
#         end
#         # println((size, count, rate))
#     end
#     approx = mean(rates[end-5:end])
#     exact = 1.5236270
#     @show (approx, exact)
#     abs(approx - exact) < 1E-2
# end

max_stars_repo_name: MiroK/fractional-dim


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Initialize data structures in element container
function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendreBasis)
  @unpack node_coordinates, left_neighbors,
          jacobian_matrix, contravariant_vectors, inverse_jacobian = elements

  linear_indices = LinearIndices(size(mesh))

  # Calculate node coordinates, Jacobian matrix, and inverse Jacobian determinant
  for cell_z in 1:size(mesh, 3), cell_y in 1:size(mesh, 2), cell_x in 1:size(mesh, 1)
    element = linear_indices[cell_x, cell_y, cell_z]

    calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, cell_z, mesh.mapping, mesh, basis)

    calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis)

    calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, node_coordinates, basis)

    calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis)
  end

  initialize_left_neighbor_connectivity!(left_neighbors, mesh, linear_indices)

  return nothing
end

# Calculate physical coordinates to which every node of the reference element is mapped
# `mesh.mapping` is passed as an additional argument for type stability (function barrier)
function calc_node_coordinates!(node_coordinates, element,
                                cell_x, cell_y, cell_z,
                                mapping, mesh::StructuredMesh{3},
                                basis::LobattoLegendreBasis)
  @unpack nodes = basis

  # Get cell length in reference mesh
  dx = 2 / size(mesh, 1)
  dy = 2 / size(mesh, 2)
  dz = 2 / size(mesh, 3)

  # Calculate node coordinates of reference mesh
  cell_x_offset = -1 + (cell_x-1) * dx + dx/2
  cell_y_offset = -1 + (cell_y-1) * dy + dy/2
  cell_z_offset = -1 + (cell_z-1) * dz + dz/2

  for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
    # node_coordinates are the mapped reference node_coordinates
    node_coordinates[:, i, j, k, element] .= mapping(cell_x_offset + dx/2 * nodes[i],
                                                     cell_y_offset + dy/2 * nodes[j],
                                                     cell_z_offset + dz/2 * nodes[k])
  end
end

# Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain
function calc_jacobian_matrix!(jacobian_matrix::AbstractArray{<:Any,6}, element, node_coordinates, basis)
  # The code below is equivalent to the following matrix multiplications but much faster.
  #
  # for dim in 1:3, j in eachnode(basis), i in eachnode(basis)
  #   # ∂/∂ξ
  #   jacobian_matrix[dim, 1, :, i, j, element] = basis.derivative_matrix * node_coordinates[dim, :, i, j, element]
  #   # ∂/∂η
  #   jacobian_matrix[dim, 2, i, :, j, element] = basis.derivative_matrix * node_coordinates[dim, i, :, j, element]
  #   # ∂/∂ζ
  #   jacobian_matrix[dim, 3, i, j, :, element] = basis.derivative_matrix * node_coordinates[dim, i, j, :, element]
  # end

  @turbo for dim in 1:3, k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
    result = zero(eltype(jacobian_matrix))

    for ii in eachnode(basis)
      result += basis.derivative_matrix[i, ii] * node_coordinates[dim, ii, j, k, element]
    end

    jacobian_matrix[dim, 1, i, j, k, element] = result
  end

  @turbo for dim in 1:3, k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
    result = zero(eltype(jacobian_matrix))

    for ii in eachnode(basis)
      result += basis.derivative_matrix[j, ii] * node_coordinates[dim, i, ii, k, element]
    end

    jacobian_matrix[dim, 2, i, j, k, element] = result
  end

  @turbo for dim in 1:3, k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
    result = zero(eltype(jacobian_matrix))

    for ii in eachnode(basis)
      result += basis.derivative_matrix[k, ii] * node_coordinates[dim, i, j, ii, element]
    end

    jacobian_matrix[dim, 3, i, j, k, element] = result
  end

  return jacobian_matrix
end

# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping,
# using the invariant curl form.
# These are called Ja^i in Kopriva's blue book.
function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any,6}, element,
                                     jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis)
  @unpack derivative_matrix = basis

  # The general form is
  # Jaⁱₙ = 0.5 * ( ∇ × (Xₘ ∇ Xₗ - Xₗ ∇ Xₘ) )ᵢ  where (n, m, l) cyclic and ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ

  for n in 1:3
    # (n, m, l) cyclic
    m = (n % 3) + 1
    l = ((n + 1) % 3) + 1

    # Calculate Ja¹ₙ = 0.5 * [ (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η - (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ ]
    # For each of these, the first and second summand are computed in separate loops
    # for performance reasons.

    # First summand 0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to j-dimension to differentiate wrt η
        result += 0.5 * derivative_matrix[j, ii] * (
          node_coordinates[m, i, ii, k, element] * jacobian_matrix[l, 3, i, ii, k, element] -
          node_coordinates[l, i, ii, k, element] * jacobian_matrix[m, 3, i, ii, k, element])
      end

      contravariant_vectors[n, 1, i, j, k, element] = result
    end

    # Second summand -0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to k-dimension to differentiate wrt ζ
        result += 0.5 * derivative_matrix[k, ii] * (
          node_coordinates[m, i, j, ii, element] * jacobian_matrix[l, 2, i, j, ii, element] -
          node_coordinates[l, i, j, ii, element] * jacobian_matrix[m, 2, i, j, ii, element])
      end

      contravariant_vectors[n, 1, i, j, k, element] -= result
    end

    # Calculate Ja²ₙ = 0.5 * [ (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ - (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ ]

    # First summand 0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to k-dimension to differentiate wrt ζ
        result += 0.5 * derivative_matrix[k, ii] * (
          node_coordinates[m, i, j, ii, element] * jacobian_matrix[l, 1, i, j, ii, element] -
          node_coordinates[l, i, j, ii, element] * jacobian_matrix[m, 1, i, j, ii, element])
      end

      contravariant_vectors[n, 2, i, j, k, element] = result
    end

    # Second summand -0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to i-dimension to differentiate wrt ξ
        result += 0.5 * derivative_matrix[i, ii] * (
          node_coordinates[m, ii, j, k, element] * jacobian_matrix[l, 3, ii, j, k, element] -
          node_coordinates[l, ii, j, k, element] * jacobian_matrix[m, 3, ii, j, k, element])
      end

      contravariant_vectors[n, 2, i, j, k, element] -= result
    end

    # Calculate Ja³ₙ = 0.5 * [ (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ - (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η ]

    # First summand 0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to i-dimension to differentiate wrt ξ
        result += 0.5 * derivative_matrix[i, ii] * (
          node_coordinates[m, ii, j, k, element] * jacobian_matrix[l, 2, ii, j, k, element] -
          node_coordinates[l, ii, j, k, element] * jacobian_matrix[m, 2, ii, j, k, element])
      end

      contravariant_vectors[n, 3, i, j, k, element] = result
    end

    # Second summand -0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η
    @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
      result = zero(eltype(contravariant_vectors))

      for ii in eachnode(basis)
        # Multiply derivative_matrix to j-dimension to differentiate wrt η
        result += 0.5 * derivative_matrix[j, ii] * (
          node_coordinates[m, i, ii, k, element] * jacobian_matrix[l, 1, i, ii, k, element] -
          node_coordinates[l, i, ii, k, element] * jacobian_matrix[m, 1, i, ii, k, element])
      end

      contravariant_vectors[n, 3, i, j, k, element] -= result
    end
  end

  return contravariant_vectors
end

# Calculate inverse Jacobian (determinant of Jacobian matrix of the mapping) in each node
function calc_inverse_jacobian!(inverse_jacobian::AbstractArray{<:Any, 4}, element, jacobian_matrix, basis)
  @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis)
    # Calculate Determinant by using Sarrus formula (about 100 times faster than LinearAlgebra.det())
    inverse_jacobian[i, j, k, element] = inv(
        jacobian_matrix[1, 1, i, j, k, element] * jacobian_matrix[2, 2, i, j, k, element] * jacobian_matrix[3, 3, i, j, k, element] +
        jacobian_matrix[1, 2, i, j, k, element] * jacobian_matrix[2, 3, i, j, k, element] * jacobian_matrix[3, 1, i, j, k, element] +
        jacobian_matrix[1, 3, i, j, k, element] * jacobian_matrix[2, 1, i, j, k, element] * jacobian_matrix[3, 2, i, j, k, element] -
        jacobian_matrix[3, 1, i, j, k, element] * jacobian_matrix[2, 2, i, j, k, element] * jacobian_matrix[1, 3, i, j, k, element] -
        jacobian_matrix[3, 2, i, j, k, element] * jacobian_matrix[2, 3, i, j, k, element] * jacobian_matrix[1, 1, i, j, k, element] -
        jacobian_matrix[3, 3, i, j, k, element] * jacobian_matrix[2, 1, i, j, k, element] * jacobian_matrix[1, 2, i, j, k, element] )
  end

  return inverse_jacobian
end

# Save id of left neighbor of every element
function initialize_left_neighbor_connectivity!(left_neighbors, mesh::StructuredMesh{3}, linear_indices)
  # Neighbors in x-direction
  for cell_z in 1:size(mesh, 3), cell_y in 1:size(mesh, 2)
    # Inner elements
    for cell_x in 2:size(mesh, 1)
      element = linear_indices[cell_x, cell_y, cell_z]
      left_neighbors[1, element] = linear_indices[cell_x - 1, cell_y, cell_z]
    end

    if isperiodic(mesh, 1)
      # Periodic boundary
      left_neighbors[1, linear_indices[1, cell_y, cell_z]] = linear_indices[end, cell_y, cell_z]
    else
      left_neighbors[1, linear_indices[1, cell_y, cell_z]] = 0
    end
  end

  # Neighbors in y-direction
  for cell_z in 1:size(mesh, 3), cell_x in 1:size(mesh, 1)
    # Inner elements
    for cell_y in 2:size(mesh, 2)
      element = linear_indices[cell_x, cell_y, cell_z]
      left_neighbors[2, element] = linear_indices[cell_x, cell_y - 1, cell_z]
    end

    if isperiodic(mesh, 2)
      # Periodic boundary
      left_neighbors[2, linear_indices[cell_x, 1, cell_z]] = linear_indices[cell_x, end, cell_z]
    else
      left_neighbors[2, linear_indices[cell_x, 1, cell_z]] = 0
    end
  end

  # Neighbors in z-direction
  for cell_y in 1:size(mesh, 2), cell_x in 1:size(mesh, 1)
    # Inner elements
    for cell_z in 2:size(mesh, 3)
      element = linear_indices[cell_x, cell_y, cell_z]
      left_neighbors[3, element] = linear_indices[cell_x, cell_y, cell_z - 1]
    end

    if isperiodic(mesh, 3)
      # Periodic boundary
      left_neighbors[3, linear_indices[cell_x, cell_y, 1]] = linear_indices[cell_x, cell_y, end]
    else
      left_neighbors[3, linear_indices[cell_x, cell_y, 1]] = 0
    end
  end

  return left_neighbors
end

end # @muladd


text:
#=
cbf:
- Julia version: 1.3.1
- Author: Berk
- Date: 2020-10-08
=#

function CBF_to_ModelData(filename; epsilon=1e-20)
    """ Converts CBF model into a fnmodel format.
        Using some code here from MathProgBase loadproblem! (technically deprecated). """
    dat = readcbfdata(filename);
    c, A, b, constr_cones, var_cones, vartypes, sense, objoffset = cbftompb(dat);
    # In MathProgBase format:
    # Conic form
    #  min  c'x
    #  st b-Ax ∈ K_1
    #        x ∈ K_2
    # Note: The sense in MathProgBase form is always minimization, and the objective offset is zero.
    # If sense == :Max, you should flip the sign of c before handing off to a solver.
    # If a cone is anything other than [:Free,:Zero,:NonNeg,:NonPos,:SOC,:SOCRotated], MPB
    # would give up. But we should be able to handle it.
    bad_cones = [:SDP];
    for (cone,idxs) in var_cones
        cone in bad_cones && error("Cone type $(cone) not supported")
    end
    n_vars = length(c);
    md = ModelData(c=c);
    # Variable bounds
    l = fill(-Inf, n_vars);
    u = fill(Inf, n_vars);
    for (cone,idxs) in var_cones
        if cone == :ExpPrimal
            l[idxs[2]] = 0;
        elseif cone == :ExpDual
            l[idxs[3]] = 0;
        elseif cone == :SOC
            l[idxs[1]] = 0;
        elseif cone == :SOCRotated
            l[idxs[1]] = 0;
            l[idxs[2]] = 0;
        else
            cone_l = (cone == :Free || cone == :NonPos) ? -Inf : 0.0;
            cone_u = (cone == :Free || cone == :NonNeg) ?  Inf : 0.0;
            for idx in idxs
                l[idx] = cone_l;
                u[idx] = cone_u;
            end
        end
    end
    # Constraint functions
    for (cone, idxs) in constr_cones
    # Idxs describe which rows of A the cone applies to.
    # All cones: [:Free, :Zero, :NonNeg, :NonPos, :SOC, :SOCRotated, :SDP, :ExpPrimal, :ExpDual]
        var_idxs = unique(cart_ind[2] for cart_ind in findall(!iszero, A[idxs, :])); # CartesianIndices...
        vks = Dict(:x => [var_idxs])
        if cone == :NonNeg
            if length(var_idxs) == 1.
                u[var_idxs[1]] = minimum([u[var_idxs[1]], b[idxs][1]/A[idxs, :].nzval[1]]);
            end
            add_lin_constr!(md, @build_constraint(A[idxs,:] * [md.vars[vk] for vk in md.vks] .<= b[idxs]));
        elseif cone == :NonPos
            if length(var_idxs) == 1.
                l[var_idxs[1]] = maximum([l[var_idxs[1]], b[idxs][1]/A[idxs, :].nzval[1]]);
            end
            add_lin_constr!(md, @build_constraint(A[idxs,:] * [md.vars[vk] for vk in md.vks] .>= b[idxs]));
        elseif cone == :Zero
            add_lin_constr!(md, @build_constraint(A[idxs,:] * [md.vars[vk] for vk in md.vks] .== b[idxs]));
        elseif cone == :SOC
            constr_fn = let b = b[idxs], A = A[idxs, var_idxs], vks = vks
                function (x)
                    vars = [x[vk] for vk in vks];
                    expr = b - A*vars;
                    return expr[1].^2 - sum(expr[2:end].^2);
                end
            end
            add_fn!(md, BlackBoxFunction(fn = constr_fn, vks = vks));
            l[var_idxs[1]] = maximum([l[var_idxs[1]], 0]);
        elseif cone == :SOCRotated
            constr_fn = let b = b[idxs], A = A[idxs, var_idxs], vks = vks
                function (x)
                    vars = [x[vk] for vk in vks];
                    expr = b - A*vars;
                    return expr[1]*expr[2] - sum(expr[3:end].^2)
                end
            end
            add_fn!(md, BlackBoxFunction(fn = constr_fn, vks = vks));
            l[var_idxs[1]] = maximum([l[var_idxs[1]], 0]);
            l[var_idxs[2]] = maximum([l[var_idxs[2]], 0]);
        elseif cone == :ExpPrimal
            constr_fn = let b = b[idxs], A = A[idxs, var_idxs], vks = vks
                function (x)
                    vars = [x[vk] for vk in vks];
                    (x,y,z) = b - A*vars;
                    return z - y*exp(x/y)
                end
            end
            add_fn!(md, BlackBoxFunction(fn = constr_fn, vks = vks));
            l[var_idxs[2]] = maximum([l[var_idxs[2]], epsilon]);
        elseif cone == :ExpDual
            constr_fn = let b = b[idxs], A = A[idxs, var_idxs], vks = vks
                function (x)
                    vars = [x[vk] for vk in vks];
                    (u,v,w) = b - A*vars;
                    if v >= 0 && w >= 0
                        return 1 # feasible
                    elseif u < 0 && w >= 0
                        return u*log(-u/w) - u + v
                    else
                        return -1 # infeasible
                    end
                end
            end
            add_fn!(md, BlackBoxFunction(fn = constr_fn, vks = vks));
            u[var_idxs[1]] = minimum([u[var_idxs[1]], 0]);
            l[var_idxs[3]] = maximum([l[var_idxs[3]], 0]);
        elseif cone in [:SDP]
            throw(ArgumentError("Haven't coded feasibility for these cones yet."));
        elseif cone == :Free
            pass;
        else
            print("This cone is not recognized.");
        end
    end
    update_bounds!(md, lbs = Dict(md.vks .=> l),
                       ubs = Dict(md.vks .=> u));
    int_vks = [Symbol("x", i) for i in findall(x -> x .== :Int, vartypes)]; # Integer var labels
    for vk in int_vks
        set_integer(md.vars[vk])
    end
    return md
end


text:
using CellListMap
using StaticArrays
import Random

#
# Here we simulate the computation of the histogram of pairwise velocities, for a set of
# particles for which the distances are known. This histogram is useful in the study
# galaxy motions, for example.
#
function pairwise_velocities(::Type{T}=Float64;
    N=100_000,cd=true,parallel=true,lcell=1,nbatches=(0,0)
) where T

    @inline dot(x::SVector{3,T}, y::SVector{3,T}) where T = x[1] * y[1] + x[2] * y[2] + x[3] * y[3]

    function compute_pairwise_mean_cell_lists!(x, y, i, j, d2, hist, velocities, rbins)
        d = x - y
        r = sqrt(d2)
        ibin = searchsortedfirst(rbins, r) - 1
        hist[1][ibin] += 1
        hist[2][ibin] += dot(velocities[i] - velocities[j], d) / r
        return hist
    end

    function reduce_hist(hist, hist_threaded)
        for i in 1:length(hist_threaded)
            hist[1] .+= hist_threaded[i][1]
            hist[2] .+= hist_threaded[i][2]
        end
        return hist
    end

    n_halos = N

    if cd
        # density galaxies in the universe megaparsec^-3
        density = T(10^5 / 20.274^3)
        boxsize = T((n_halos/density)^(1/3))
    else
        boxsize = T(20.274)
    end

    Random.seed!(321)
    Lbox = T[boxsize,boxsize,boxsize]
    positions = convert.(T,boxsize .* rand(Float64, 3, n_halos))
    velocities = convert.(T,rand(Float64, 3, n_halos))
    rbins = T[0.,1.,2.,3.,4.,5.]
    r_max = maximum(rbins)

    n = size(positions)[2]
    positions = reshape(reinterpret(SVector{3,T}, positions), n)
    velocities = reshape(reinterpret(SVector{3,T}, velocities), n)

    box = Box(Lbox, r_max, UnitCellType=OrthorhombicCell, lcell=lcell)
    cl = CellList(positions, box, parallel=parallel, nbatches=nbatches)
    hist = (zeros(Int, length(rbins) - 1), zeros(T, length(rbins) - 1))

    # Needs this to stabilize the type of velocities and hist, probably
    function barrier!(f!, velocities, rbins, hist, box, cl, reduce_hist, parallel)
        map_pairwise!(
            (x, y, i, j, d2, hist) -> f!(x, y, i, j, d2, hist, velocities, rbins),
            hist, box, cl,
            reduce=reduce_hist,
            parallel=parallel
        )
        return hist
    end

    barrier!(
        compute_pairwise_mean_cell_lists!,
        velocities,rbins,hist,box,cl,reduce_hist,parallel
    )

    n_pairs = hist[1]
    mean_v_r = hist[2]
    mean_v_r[n_pairs .> 0] = mean_v_r[n_pairs .> 0] ./ n_pairs[n_pairs .> 0]

    # using Julia 1.7.0, this is the expected result
    correct_result = [
        0.00028420429835911875,
       -3.240572326959191e-5,
        9.581744218539093e-5,
       -2.879083693540441e-6,
       -3.8499932402028564e-5
       ]
    return mean_v_r ≈ correct_result, mean_v_r

end

# Credits to @florpi, and we use her github name here as an alias for testing
# and historical reasons.
const florpi = pairwise_velocities


text:
function hmc_step(current_state::state, log_likelihood::Function, grad_log_likelihood::Function, current_p::Vector{<:Real}, ϵ::Float64, L::Int)

    current_log_likelihood = current_state.log_likelihood[1] + 0.5 * sum(current_p.^2)

    # Evolve p half-step
    current_p += 0.5 * ϵ * current_state.grad_log_likelihood

    # Leap-Frog Integrate
    x1 = current_state.x
    for i = 1:L-1
        x1 += ϵ * current_p
        current_p += ϵ * grad_log_likelihood(x1)
    end
    x1 += ϵ * current_p

    current_p += 0.5 * ϵ * grad_log_likelihood(x1)
    current_p = -1 * current_p

    # Perform comparison
    new_log_likelihood = log_likelihood(x1)
    proposed_log_likelihood = new_log_likelihood + 0.5 * sum(current_p.^2)
    alpha = current_log_likelihood - proposed_log_likelihood
    u = log(rand())
    if u <= alpha
        current_state.x[:] .= x1
        current_state.log_likelihood[1] = new_log_likelihood
        return 1
    else
        return 0
    end

end

function hmc(x0 :: Vector{<:Real}, log_likelihood :: Function, grad_likelihood_function :: Function, num_burnin :: Int, num_samples :: Int)

    # Initialize
    ϵ = 1e-2
    L = 20

    # Initialize
    current_state = state(x0, log_likelihood(x0), grad_likelihood_function(x0))

    # Burnin
    burnin_chain = chain(x0, num_burnin)
    burnin_chain.diagnostics.log_likelihood[1] = log_likelihood(x0)
    for iteration = 1:num_burnin-1
        p = randn(size(x0))
        flag = hmc_step(current_state, log_likelihood, grad_likelihood_function, p, ϵ, L)
        burnin_chain.x[iteration+1] = copy(current_state.x)
        burnin_chain.diagnostics.log_likelihood[iteration+1] = current_state.log_likelihood[1]
        burnin_chain.diagnostics.acceptance_ratio[iteration+1] = ( burnin_chain.diagnostics.acceptance_ratio[iteration] * (iteration - 1) + flag ) / iteration
    end

    if burnin_chain.diagnostics.acceptance_ratio[end] < 0.25
        ϵ *= 0.5
    elseif burnin_chain.diagnostics.acceptance_ratio[end] > 0.9
        ϵ *= 2
    end

    # Sampling
    sampling_chain = chain(current_state.x, num_samples)
    sampling_chain.diagnostics.log_likelihood[1] = log_likelihood(x0)

    for iteration = 1:num_samples-1
        p = randn(size(x0))
        flag = hmc_step(current_state, log_likelihood, grad_likelihood_function, p, ϵ, L)
        sampling_chain.x[iteration+1] = copy(current_state.x)
        sampling_chain.diagnostics.log_likelihood[iteration+1] = current_state.log_likelihood[1]
        sampling_chain.diagnostics.acceptance_ratio[iteration+1] = ( sampling_chain.diagnostics.acceptance_ratio[iteration] * (iteration - 1) + flag ) / iteration
    end

    return burnin_chain, sampling_chain
end


text:
module TestAqua

import Aqua
import SyncBarriers

test() = Aqua.test_all(SyncBarriers; unbound_args = false)

end  # module


text:
module BondGraphs

import Graphs as g
import Base: RefValue, eltype, show, in, iterate, ==, getproperty
import ModelingToolkit: parameters, states, equations # Importing means names can be reused, but may be confusing

using StaticArrays
using ModelingToolkit
using DifferentialEquations
using SymbolicUtils, SymbolicUtils.Rewriters
using OrderedCollections
using Setfield
using Catalyst
using RecipesBase, GraphRecipes

export AbstractNode, Component, Junction, EqualEffort, EqualFlow,
SourceSensor, Port, Bond, BondGraph, BondGraphNode,

type, name, freeports, numports, weights, vertex, set_vertex!,
parameters, states, defaults, constitutive_relations,
get_parameter, set_parameter!, get_initial_value, set_initial_value!,

srcnode, dstnode, nodes, bonds, components, junctions, getnodes, getbonds,

add_node!, remove_node!, connect!, disconnect!,
swap!, insert_node!, merge_nodes!, simplify_junctions!, expose,

simulate, addlibrary!

# Component libraries
include("libraries/biochemical.jl")
include("libraries/standard.jl")
include("libraries/libraryfunctions.jl")

# Types used by BondGraphs
include("basetypes/AbstractNode.jl")
include("basetypes/Bond.jl")
include("basetypes/BondGraph.jl")

# Core functionality
include("graphfunctions.jl")
include("construction.jl")
include("equations.jl")
include("conversion.jl")
include("plotrecipes.jl")

end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
@testset "Dynamic Oracles" begin

    oracle = Oracle(ArcEager(), dynamic_oracle, typed)

    model(cfg) = nothing

    graph = test_sentence("economicnews.conll")

    model(cfg) = static_oracle(cfg, graph, typed)

    cfg = initconfig(oracle, graph)
    nocost(t, cfg) = DependencyTrees.TransitionParsing.cost(t, cfg, graph) == 0
    while !isfinal(cfg)
        pred = model(cfg)
        G = oracle(cfg, graph)
        @test pred in G
        @test all(t -> nocost(t, cfg), G)
        cfg = pred(cfg)
    end

    # make sure this works the same for untyped oracles too
    oracle = Oracle(ArcEager(), dynamic_oracle, untyped)
    model(cfg) = static_oracle(cfg, graph, untyped)
    cfg = initconfig(oracle.system, graph)
    while !isfinal(cfg)
        pred = model(cfg)
        G = oracle(cfg, graph)
        @test pred in G
        @test any(t -> !nocost(t, cfg), [Shift(), Reduce(), LeftArc(), RightArc()])
        @test all(t -> nocost(t, cfg), G)
        cfg = pred(cfg)
    end

    tbfile = joinpath(@__DIR__, "data", "wsj_0001.dp")
    readtok = line -> begin
        form, deprel, head = split(line, "\t")
        head = parse(Int, head)
        DependencyTrees.Token(form, head, deprel)
    end
    treebank = Treebank(tbfile, readtok)

    for (cfg, gold) in oracle(graph)
        @test length(gold) >= 1
    end

    @testset "Projectivity" begin
        o = Oracle(ArcHybrid(), dynamic_oracle)
        tb = Treebank(joinpath(@__DIR__, "data", "nonprojective1.conll"))
        @test length(collect(tb)) == 1
        @test length(collect(Iterators.flatten(oracle.(tb)))) == 0
        @test sum(length(o(t)) for t in tb) == 0

        for tree in tb
            if !is_projective(tree)
                @test isempty(oracle(tree))
            else
                @test !isempty(oracle(tree))
            end
        end
    end

    @testset "Dynamic Iteration" begin
        oracle = Oracle(ArcHybrid(), dynamic_oracle, untyped)
        policy = NeverExplore()

        tb = Treebank(joinpath(@__DIR__, "data", "hybridtests.conll"))

        for tree in treebank
            for (cfg, G) in oracle(tree)
                @test cfg isa DependencyTrees.TransitionParsing.ArcHybridConfig
                for t in G
                    T = typeof(t)
                    @test T <: LeftArc || T <: RightArc || T <: Shift
                end
            end

            for (cfg, G) in oracle(tree)
                A = possible_transitions(cfg, tree)
                @test G ⊆ A
                t = policy(cfg, A, G)
                @test t in A
            end

            # explore transition space manually
            gold = oracle(tree)
            cfg = initconfig(oracle, tree)
            while !isfinal(cfg)
                t = TransitionParsing.choose_transition(gold, cfg)
                A = TransitionParsing.possible_transitions(cfg)
                @test t in A
                cfg = t(cfg)
            end
        end
    end
end

@testset "Transition Output Spaces" begin
    @test transition_space(ArcEager()) == [LeftArc(), RightArc(), Reduce(), Shift()]
    @test transition_space(ArcEager(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                      RightArc("a"), RightArc("b"),
                                                      Reduce(), Shift()]

    @test transition_space(ArcHybrid()) == [LeftArc(), RightArc(), Shift()]
    @test transition_space(ArcHybrid(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                       RightArc("a"), RightArc("b"),
                                                       Shift()]

    @test transition_space(ArcStandard()) == [LeftArc(), RightArc(), Shift()]
    @test transition_space(ArcStandard(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                         RightArc("a"), RightArc("b"),
                                                         Shift()]

    @test transition_space(ArcSwift(); max_k=2) == [LeftArc(1),  LeftArc(2),
                                                    RightArc(1), RightArc(2),
                                                    Shift()]
    ts1 = Set(transition_space(ArcSwift(), ["a","b"]; max_k=2))
    ts2 = Set([LeftArc(1, "a"),  LeftArc(2, "a"),
               LeftArc(1, "b"),  LeftArc(2, "b"),
               RightArc(1, "a"), RightArc(2, "a"),
               RightArc(1, "b"), RightArc(2, "b"),
               Shift()])
    @test  ts1 == ts2

    @test transition_space(ListBasedNonProjective()) == [LeftArc(), RightArc(), NoArc(), Shift()]
    @test transition_space(ListBasedNonProjective(), ["a","b"]) == [LeftArc("a"), LeftArc("b"),
                                                                    RightArc("a"), RightArc("b"),
                                                                    NoArc(), Shift()]
end

@testset "Exploration Policies" begin
    always1, never1 = AlwaysExplore(), NeverExplore()
    always2, never2 = ExplorationPolicy(1), ExplorationPolicy(0)
    @test showstr(always1) == "AlwaysExplore"
    @test showstr(never1) == "NeverExplore"
    for i = 1:10
        @test always1() == always2() == true
        @test never1()  == never2()  == false
    end

    @testset "Models" begin
        oracle = Oracle(ArcEager(), dynamic_oracle, untyped)
        tree = first(test_treebank("chopsticks.conll"))
        @testset "Always" begin
            policy = AlwaysExplore((cfg, A, G) -> first(A))
            xys = collect(oracle(tree, policy))
        end
        @testset "Never" begin
            policy = NeverExplore((cfg, A, G) -> first(A))
            xys = collect(oracle(tree, policy))
        end
        @testset "Rate" begin
            policy = ExplorationPolicy(0.1, (cfg, A, G) -> first(A))
            xys = collect(oracle(tree, policy))
        end
    end

    @testset "RNG" begin
        using Random
        rng = MersenneTwister(42)
        always = AlwaysExplore(rng)
        never  = NeverExplore(rng)
        sometimes = ExplorationPolicy(0.5, rng)

        tree = test_sentence("economicnews.conll")
        for system in (ArcEager(), ArcHybrid())
            oracle = Oracle(system, dynamic_oracle)
            for policy in (always, never, sometimes)
                for (cfg, G) in oracle(tree, policy)
                    @test length(G) >= 1
                end
            end
        end
    end

end

@testset "Feature Extraction" begin

    function check(oracle, gold_tree)
        for (cfg, G) in oracle(gold_tree)
            stk, buf = stack(cfg), buffer(cfg)
            ts = [buffertoken(cfg, b) for b in buf]
            ts = [stacktoken(cfg, s) for s in stk]
        end
        return true
    end

    treebank = Treebank(joinpath(@__DIR__, "data", "hybridtests.conll"))
    for system in (ArcEager(), ArcHybrid())
        for oracle_fn in (static_oracle, dynamic_oracle)
            oracle = Oracle(system, oracle_fn, untyped)
            @test all(check(oracle, tree) for tree in treebank)
        end
    end
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSLogFull

using Test, RLinearAlgebra, Random, LinearAlgebra

Random.seed!(1010)

@testset "LS Full Log -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LSLogFull) == LinSysSolverLog

    # Verify Required Fields
    @test :iterations in fieldnames(LSLogFull)
    @test :converged in fieldnames(LSLogFull)

    # Verify log_update initialization
    let
        A = rand(2,2)
        x = rand(2)
        b = A * x
        z = rand(2)

        sampler = LinSysVecRowOneRandCyclic()
        log = LSLogFull()

        RLinearAlgebra.log_update!(log, sampler, z, (), 0, A, b)

        @test length(log.resid_hist) == 1
        @test log.resid_hist[1] == norm(A * z - b)
        @test log.iterations == 0
        @test log.converged == false
    end

    # Verify collection rate
    let
        A = rand(2,2)
        x = rand(2)
        b = A * x
        z = rand(2)

        sampler = LinSysVecRowOneRandCyclic()
        log = LSLogFull(3)

        # Test initialization of log
        RLinearAlgebra.log_update!(log, sampler, z, (), 0, A, b)

        @test length(log.resid_hist) == 1
        @test log.resid_hist[1] == norm(A * z - b)
        @test log.iterations == 0
        @test log.converged == false

        # Test collection rate of log
        for i = 1:10
            RLinearAlgebra.log_update!(log, sampler, x + i*(z-x), (), i, A, b)
        end

        @test length(log.resid_hist) == 4 # Record at 0, 3, 6, 9
        @test norm(log.resid_hist - [1.0, 3.0, 6.0, 9.0] * norm(A * z - b)) < 1e-15
        @test log.iterations == 10
        @test log.converged == false
    end

end

end # End Module


text:
#=
Small illustation on the capacitated vehicle routing problem (CVRP)
Just run the script for an illustration

Mathematical model based on: https://www.researchgate.net/publication/323173028_Two_models_of_the_capacitated_vehicle_routing_problem
using modified indices to match Julia's 1-based indexing

In this model we eliminate all subtour explicitely, keep in mind that this is very expensive
as the number of subtours grows exponentially with the number of cities! A modified approach is to solve the
relaxed problem (i.e. without subtour constraints) and only add them when they occur (thus solving the model again),
adding the constraint that prohibits the subtour that was found during the previous iteration cf.
https://how-to.aimms.com/Articles/332/332-Implicit-Dantzig-Fulkerson-Johnson.html
for a clear explanation.

=#
using JuMP
using GLPK
using Plots
using Combinatorics:powerset # may require additional install

"""
Multiple Vehicle Routing problem (using the entire set of possible subtours for subtour elimination)
"""
function MVRP()
    # general problem setup #
    # --------------------- #
    # number of cities
    n = 6
    # pick random locations
    coord_x = rand(n+1)
    coord_y = rand(n+1)
    dist(x_1, x_2, y_1, y_2) = sqrt((x_1 - x_2)^2 + (y_1 - y_2)^2)
    # cost matrix (symmetrical, distances, depot is city 1)
    c = zeros(n+1, n+1)
    for i in 1:n+1
        for j in i+1:n+1
            c[i,j] = dist(coord_x[i], coord_x[j], coord_y[i], coord_y[j])
        end
    end
    c .+= c'
    # demand for each city (number of e.g. number of pallets)
    d = rand(1:10, n+1)
    d[1] = 0 # no demand @ depot
    # truck capacity (equal for all trucks, can easily be extended to different size trucks)
    Q = 20
    # minimum number of trucks for order to be feasible:
    p = ceil(Int,sum(d) / Q)
    @info "$(p) trucks required"

    # Linear programming setup #
    # ------------------------ #
    model = Model(with_optimizer(GLPK.Optimizer))
    # Add variables (x_{i,j,p}: edge from i > j by truck p, Boolean valued)
    @variable(model, x[1:n+1,1:n+1,1:p], Bin)
    # Add objective
    @objective(model, Min, sum(sum(x[:,:,i] .* c) for i in 1:p))
    # Add constraints
    for r in 1:p
        # each vehicle can leave the depot only once
        @constraint(model, sum(x[1, j, p] for j in 2:n+1) == 1)
        # sum of the demands of the customers visited in a route is less than or equal to the capacity of the vehicle
        @constraint(model, sum(d[j]*x[i,j,r] for i in 1:n+1 for j in 2:n+1 if i≠j) <= Q)
        # the number of the vehicles arriving at every customer and entering the depot is equal to the number of the vehicles leaving
        for j in 1:n+1
            @constraint(model, sum(x[i,j,r] for i in 1:n+1 if i≠j) == sum(x[j,i,r] for i in 1:n+1))
        end
    end
    # ensure that each customer is visited by exactly one vehicle
    for j in 2:n+1
        @constraint(model, sum(x[i,j,r] for r in 1:p for i in 1:n+1 if i≠j) == 1)
    end

    #=
    # eliminate sub-tours (essential but costly)
    for subset in powerset(2:n+1, 2, n)
        @constraint(model, sum(x[i,j,r] for r in 1:p for i in subset for j in subset if i≠j) <= length(subset) - 1 )
    end
    =#
    # Solve problem
    optimize!(model)
    @info "termination status: $(termination_status(model))"
    @info "objective_value $(objective_value(model))"

    #     Nice illustration    #
    # ------------------------ #
    marks = [:circle for _ in 1:n+1]
    marks[1] = :cross
    fig = plot()
    scatter!(fig, coord_x[1:1], coord_y[1:1], marker=:cross, label="Depot", color=:black,markersize=12)
    scatter!(fig, coord_x[2:end], coord_y[2:end], marker=:circle, label="Stores", color=:gray )
    for r in 1:p
        inds = findall(x-> !iszero(x), value.(x[:,:,r]))
        for ind in inds
            xx = [coord_x[ind[1]]; coord_x[ind[2]]]
            yy = [coord_y[ind[1]]; coord_y[ind[2]]]
            plot!(fig, xx, yy, color=r, label="")
        end
    end
    title!(fig, "lay-out for $(p) trucks with $(n) cities")
    display(fig)
end

"""
Multiple Trip Vehicle Routing problem

Implementation according to the 4-index formulation, also uses the entire set of subtours (similar to the other example)

Remark: since we are mimimising the total time of all trips, it can happen that a single truck gets two trips, while the
other one gets nothing at all. In a simulation context, you should take this into account.
"""
function MTVRP()
    # general problem setup #
    # --------------------- #
    # number of cities
    n = 7
    # pick random locations
    coord_x = rand(n+1)
    coord_y = rand(n+1)
    dist(x_1, x_2, y_1, y_2) = sqrt((x_1 - x_2)^2 + (y_1 - y_2)^2)
    # cost matrix (symmetrical, distances, depot is city 1)
    c = zeros(n+1, n+1)
    for i in 1:n+1
        for j in i+1:n+1
            c[i,j] = dist(coord_x[i], coord_x[j], coord_y[i], coord_y[j])
        end
    end
    c .+= c'
    # demand for each city (number of e.g. number of pallets)
    d = rand(1:10, n+1)
    d[1] = 0 # no demand @ depot
    # truck capacity (equal for all trucks, can easily be extended to different size trucks)
    Q = 20
    V = 2 # number of vehicles
    R = 3 # number of trips/vehicle
    TH = 5 # max trip length => results in "infeasible"
    # minimum number of trucks for order to be feasible:
    @info "$(ceil(Int,sum(d) / Q)) trips required for a solution"
    @info "$(V) trucks used, max $(R) trips per truck, total work time $(TH)"
    # Linear programming setup #
    # ------------------------ #
    model = Model(with_optimizer(GLPK.Optimizer))
    # Add variables
    # x_{i,j,v, r}: edge from i > j by truck v on trip r, Boolean valued
    @variable(model, x[1:n+1, 1:n+1, 1:V, 1:R], Bin)
    # y_{i, v, r}: 1 if trip r of vehicle v visits vertex i, Boolean valued
    @variable(model, y[1:n+1, 1:V, 1:R], Bin)

    # Add objective (slices)
    @objective(model, Min, sum(sum(x[:,:,v,r] .* c) for v in 1:V for r in 1:R)) # (1)

    # Add constraints
    for i in 2:n+1
        @constraint(model, sum(y[i,v,r] for v in 1:V for r in 1:R) == 1) # (2)
    end

    for v in 1:V
        for r in 1:R
            for i in 1:n+1
                @constraint(model, sum(x[i,j,v,r] for j in 1:n+1 if j≠i) == y[i,v,r]) # (3)
                @constraint(model, sum(x[j,i,v,r] for j in 1:n+1 if j≠i) == y[i,v,r]) # (3)
                @constraint(model, sum(x[j,i,v,r] for j in 1:n+1 if j≠i) == sum(x[i,j,v,r] for j in 1:n+1 if j≠i)) # (3)
            end

            @constraint(model, sum(d[i] * y[i,v,r] for i in 2:n+1) <= Q) # (4)
            for subset in powerset(2:n+1, 2, n)
                @constraint(model, sum(x[i,j,v,r] for i in subset for j in subset if i≠j) <= length(subset) - 1 ) # (5)
            end

        end
    end

    for v in 1:V
        @constraint(model, sum(sum(x[:,:,v,r] .* c) for r in 1:R) <= TH) # (6)
    end

    optimize!(model)
    @info "termination status: $(termination_status(model))"
    @info "objective_value $(objective_value(model))"

    #     Nice illustration    #
    # ------------------------ #
    marks = [:circle for _ in 1:n+1]
    marks[1] = :cross
    fig = plot()
    scatter!(fig, coord_x[1:1], coord_y[1:1], marker=:cross, label="Depot", color=:black,markersize=12)
    scatter!(fig, coord_x[2:end], coord_y[2:end], marker=:circle, label="Stores", color=:gray )
    for v in 1:V
        for r in 1:R
            inds = findall(x-> !iszero(x), value.(x[:,:,v,r]))
            for ind in inds
                xx = [coord_x[ind[1]]; coord_x[ind[2]]]
                yy = [coord_y[ind[1]]; coord_y[ind[2]]]
                plot!(fig, xx, yy, color=v, label="", line=r)
            end
            @info "truck $(v), route $(r): $(round(sum(c[inds]), digits=2)) kms"
        end
    end
    title!(fig, "lay-out for $(V) trucks with $(n) cities, max $(R) routes per truck\ncolor= truck, line thickness=route")
    display(fig)
end

MVRP()
#MTVRP()


text:
#=
# 405: Generic Operator: 1D Nonlinear Poisson equation
([source code](SOURCE_URL))

Solve the nonlinear Poisson equation

```math
-\nabla \varepsilon \nabla u + e^{u}-e^{-u} = f
```
in $\Omega=(0,1)$ with boundary condition $u(0)=0$ and $u(1)=1$ with
```math
f(x)=
    \begin{cases}
    1&,x>0.5\\
    -1&, x<0.5
    \end{cases}.
```

This stationary problem is an example of a nonlinear Poisson equation or Poisson-Boltzmann equation.
Such equation occur e.g. in simulations of electrochemical systems and semicondutor devices.

=#

module Example405_GenericOperator
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse)

    ## Create a one-dimensional discretization
    h=1.0/convert(Float64,n)
    X=collect(0:h:1)
    grid=VoronoiFVM.Grid(X)

    ## A parameter which is "passed" to the flux function via scope
    ϵ=1.0e-3

    ## This generic operator works on the full solution seen as linear vector, and indexing
    ## into it needs to be performed with the help of idx (defined below for a solution vector)
    ## Here, instead of the flux function we provide a "generic operator"
    ## which provides the stiffness part of the problem. Its sparsity is detected automatically
    ## using SparsityDetection.jl
    function generic_operator!(f,u,sys)
        f.=0.0
        for i=1:length(X)-1
            du=ϵ*(u[idx[1,i]]-u[idx[1,i+1]])/(X[i+1]-X[i])
            f[idx[1,i]]+=du
            f[idx[1,i+1]]-=du
        end
    end

    ## Source term
    function source!(f,node)
        if node[1]<=0.5
            f[1]=1
        else
            f[1]=-1
        end
    end

    ## Reaction term
    function reaction!(f,u,node)
        f[1]=exp(u[1]) - exp(-u[1])
    end

    ## Create a physics structure
    physics=VoronoiFVM.Physics(
        generic=generic_operator!,
        source=source!,
        reaction=reaction!)

    ## Create a finite volume system - either
    ## in the dense or  the sparse version.
    ## The difference is in the way the solution object
    ## is stored - as dense or as sparse matrix

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    ## Add species 1 to region 1
    enable_species!(sys,1,[1])

    ## Set boundary conditions
    boundary_dirichlet!(sys,1,1,0.0)
    boundary_dirichlet!(sys,1,2,1.0)

    ## Create a solution array
    inival=unknowns(sys,inival=0.5)
    solution=unknowns(sys)

    idx=unknown_indices(solution)
    ## Create solver control info
    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose

    ## Stationary solution of the problem
    solve!(solution,inival,sys, control=control)

    scalarplot(grid,solution[1,:],title="Nonlinear Poisson",Plotter=Plotter)
    return sum(solution)
end

function test()
    testval=1.5247901344230088
    main(unknown_storage=:sparse) ≈ testval && main(unknown_storage=:dense) ≈ testval
end

end


text:
abstract type AbstractWeights end
using DataFrames
mutable struct Weights <: AbstractWeights
    lookup::Dict
end
mutable struct NoWeights <: AbstractWeights
    lookup::Dict
end
function NormalWeight(df::DataFrame, y::Array{Any})
    features = eachcol(df)
    n_features = size(df)[1]
    categories = Set(y)
    y = OrdinalEncoder(y).predict(y)
    dist = NormalDist(y)
    norm = dist.apply(y)
    for val in y

    end
end
function ManualWeights(pairs::Array)
    return(Weights(Dict(pairs)))
end


text:
using Oceananigans.Grids: Center, Face
using Oceananigans.Fields: XFaceField, YFaceField, ZFaceField
using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions, VVelocityBoundaryConditions, FieldBoundaryConditions

function HydrostaticFreeSurfaceVelocityFields(::Nothing, arch, grid, clock, bcs=NamedTuple())
    u_bcs = :u ∈ keys(bcs) ? bcs.u : UVelocityBoundaryConditions(grid)
    v_bcs = :v ∈ keys(bcs) ? bcs.v : VVelocityBoundaryConditions(grid)
    w_bcs = :w ∈ keys(bcs) ? bcs.w : FreeSurfaceWVelocityBoundaryConditions(grid)

    u = XFaceField(arch, grid, u_bcs)
    v = YFaceField(arch, grid, v_bcs)
    w = ZFaceField(arch, grid, w_bcs)

    return (u=u, v=v, w=w)
end

FreeSurfaceWVelocityBoundaryConditions(grid; user_defined_bcs...) =
    FieldBoundaryConditions(grid, (Center, Center, Face); top=nothing, user_defined_bcs...)


text:
import DualNumbers

export ∇n

function gradient(f::Function,Γ::Domain)
    fx = Fun(x->f(DualNumbers.Dual(real(x),1),imag(x)),Γ)
    fy = Fun(x->f(real(x),DualNumbers.Dual(imag(x),1)),Γ)
    Fun(Γ,DualNumbers.epsilon.(fx.coefficients)),Fun(Γ,DualNumbers.epsilon.(fy.coefficients))
end

∇n(f::Function,d::Domain) = normalderivative(f,d)

normalderivative(f::Function,Γ::UnionDomain) = Fun(map(d->normalderivative(f,d),Γ),PiecewiseSpace)

function normalderivative(f::Function,Γ::Domain)
    fx,fy = gradient(f,Γ)
    rp = Fun(x->fromcanonical(Γ,x),Γ)'
    n = -im*rp/abs(rp)
    fx*real(n)+fy*imag(n)
end


text:
# This module overwrite the default graphic width
# for the sake of comfy-ness
using Gadfly
using Gadfly.Compose
using Cairo
using Fontconfig
using Distributions

const DEF_WIDTH=Compose.default_graphic_width
const DEF_HEIGHT=Compose.default_graphic_height

function qqplot(X; transpose=false, w=DEF_WIDTH, h=DEF_WIDTH)
	set_default_plot_size(w, h)

	d = Normal(0, 1)
	n = length(X)
	sx = sort(X)
	q = map(eachindex(X)) do i
		quantile(d, (i - 0.5) / n)
	end
	theme = Theme(
		discrete_highlight_color=x -> "teal",
		default_color="ivory"
	)
	xlabel = "𝐪"
	ylabel = "𝐫"

	if transpose
		sx, q = q, sx
		xlabel, ylabel = ylabel, xlabel
	end
	points = layer(x=q, y=sx, Geom.point, theme, )
	lines = layer(x=q, y=sx, Geom.line, color=[colorant"teal"])
	plot(lines, points, Guide.xlabel(xlabel), Guide.ylabel(ylabel))
end

function coefplot(model; w=DEF_WIDTH, bar_height = 0.5cm, intercept=false, sortabs=false, filter=true, limit=nothing, alpha=0.05)
	abs_col = :AbsCoef
	coef_col = "Coef."
	df = DataFrame(coeftable(model))
	n = size(modelmatrix(model), 1)
	t = cquantile(TDist(n-2), alpha/2)
	codf = @chain begin
		df
		# Remove 95%
		if filter
			Base.filter(r ->  !(-t ≤ r.t ≤ t), _)
			Base.filter(r -> r["Coef."] ≉ 0, _)
		else
			_
		end
		# Remove intercept term
		if intercept
			_
		else
			Base.filter(r->r["Name"] != "(Intercept)", _)
		end
		# abs value of coef column
		transform(_, coef_col => ByRow(abs) => abs_col)
		# take from top
		if isnothing(limit)
			_
		else
			df = sort(_, coef_col, rev=true)
			k = limit÷2
			[first(df, k); last(df, k)]
		end
		if sortabs
			sort(_, abs_col, rev=false)
		else
			sort(_, coef_col, rev=false)
		end
	end
	set_default_plot_size(w, bar_height * size(codf, 1))
	theme = Theme(bar_highlight=colorant"teal", default_color="teal")
	geo = Geom.bar(orientation=:horizontal)
	p = plot(codf, x="Coef.", y="Name", Guide.xlabel("Value"), Guide.ylabel("Coef"), geo, theme)
end

function plot_pred_eps(pred, e; transpose=false, w=DEF_WIDTH, h=DEF_WIDTH)
	set_default_plot_size(w, h)

	theme = Theme(
		discrete_highlight_color=x -> "teal",
		default_color="ivory"
	)
	xlabel = "𝐲"
	ylabel = "𝛜"
	x, y = pred, e

	if transpose
		x, y = y, x
		xlabel, ylabel = ylabel, xlabel
	end
	points = layer(x=x, y=y, Geom.point, theme)
	plot(points, Guide.xlabel(xlabel), Guide.ylabel(ylabel))
end

function plot_eps_hist(e; w=DEF_WIDTH, h=DEF_HEIGHT)

	set_default_plot_size(w, h)

	theme = Theme(
		discrete_highlight_color=x -> "teal",
		default_color="teal"
	)
	plot(x=e, y=e, Guide.xlabel("𝐫"), Guide.ylabel(""), Geom.histogram, theme)
end


text:
function venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, idx, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
    # @cuprintf("Type is %s", typeof(VL_CONST))
    epsigh = gpuConfigData[8] * gpuGlobalDataFixedPoint[idx].short_distance
    epsi = epsigh*epsigh*epsigh

    # shared[thread_idx], shared[thread_idx + block_dim * 1], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] = 1,1,1,1

    for i in 0:3
        shared[thread_idx + block_dim * i] = 1
        q = gpuGlobalDataRest[9+i, idx]
        del_neg = gpuGlobalDataRest[9+i, idx] - 0.5*(delx * gpuGlobalDataRest[13+i, idx] + dely * gpuGlobalDataRest[17+i, idx]) - q
        if abs(del_neg) > 1e-5
            del_pos = gpuGlobalDataRest[21+i, idx] - q
            if del_neg < 0
                del_pos = gpuGlobalDataRest[25+i, idx] - q
            end
            if del_neg == 0
                del_pos = 0
            end
            num = (del_pos*del_pos) + (epsi*epsi)
            num = (num*del_neg) + 2 * (del_neg*del_neg*del_pos)

            den = (del_pos*del_pos) + (2 *del_neg*del_neg)
            den += (del_neg*del_pos) + (epsi*epsi)
            den *= del_neg

            temp = num/den
            if temp > 1
                temp = 1
            end
            shared[thread_idx + block_dim * i] = temp
        end
    end
    return nothing
end

@inline function qtilde_to_primitive_kernel(qtilde, gpuConfigData, shared, thread_idx, block_dim)
    # # gamma = gpuConfigData[15]
    beta = -qtilde[4]*0.5
    temp = 0.5/beta
    u1 = qtilde[2]*temp
    u2 = qtilde[3]*temp

    temp2 = qtilde[1] + beta*(u1*u1 + u2*u2) - (CUDAnative.log(beta)/(gpuConfigData[15]-1))
    # rho = CUDAnative.exp(temp2)
    shared[thread_idx + block_dim * 4] = u1
    shared[thread_idx + block_dim * 5] = u2
    shared[thread_idx + block_dim * 6] = CUDAnative.exp(temp2)
    shared[thread_idx + block_dim * 7] = shared[thread_idx + block_dim * 6]*temp
    return nothing
end

# function max_q_values_kernel(gpuGlobalDataCommon, idx, maxq)
#     maxq = (
#                 gpuGlobalDataRest[9, idx],
#                 gpuGlobalDataRest[10, idx],
#                 gpuGlobalDataRest[11, idx],
#                 gpuGlobalDataRest[12, idx]
#             )
#     for iter in 9:28
#         conn = Int(gpuGlobalDataCommon[iter, idx])
#         if conn == 0.0
#             break
#         end
#         # currq = globaldata[itm].q
#         if maxq[1] < gpuGlobalDataRest[9, conn]
#             maxq = (gpuGlobalDataRest[9, conn], maxq[2], maxq[3], maxq[4])
#         end
#         if maxq[2] < gpuGlobalDataRest[10, conn]
#             maxq = (maxq[1], gpuGlobalDataRest[10, conn], maxq[3], maxq[4])
#         end
#         if maxq[3] < gpuGlobalDataRest[11, conn]
#             maxq = (maxq[1], maxq[2], gpuGlobalDataRest[11, conn], maxq[4])
#         end
#         if maxq[4] < gpuGlobalDataRest[12, conn]
#             maxq = (maxq[1], maxq[2], maxq[3], gpuGlobalDataRest[12, conn])
#         end
#     end
#     return nothing
# end

# function min_q_values_kernel(gpuGlobalDataCommon, idx, minq)
#     minq = (
#                 gpuGlobalDataRest[9, idx],
#                 gpuGlobalDataRest[10, idx],
#                 gpuGlobalDataRest[11, idx],
#                 gpuGlobalDataRest[12, idx]
#             )
#     for iter in 9:28
#         conn = Int(gpuGlobalDataCommon[iter, idx])
#         if conn == 0.0
#             break
#         end
#         # currq = globaldata[itm].q
#         if minq[1] > gpuGlobalDataRest[9, conn]
#             minq = (gpuGlobalDataRest[9, conn], minq[2], minq[3], minq[4])
#         end
#         if minq[2] > gpuGlobalDataRest[10, conn]
#             minq = (minq[1], gpuGlobalDataRest[10, conn], minq[3], minq[4])
#         end
#         if minq[3] > gpuGlobalDataRest[11, conn]
#             minq = (minq[1], minq[2], gpuGlobalDataRest[11, conn], minq[4])
#         end
#         if minq[4] > gpuGlobalDataRest[12, conn]
#             minq = (minq[1], minq[2], minq[3], gpuGlobalDataRest[12, conn])
#         end
#     end
#     return nothing
# end

# @inline function smallest_dist(globaldata, idx::Int64)
#     min_dist = 1000.0
#     for itm in globaldata[idx].conn
#         ds = hypot(globaldata[idx].x - globaldata[itm].x, globaldata[idx].y - globaldata[itm].y)
#         if ds < min_dist
#             min_dist = ds
#         end
#     end
#     globaldata[idx].short_distance = min_dist
# end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
function twodMassMatrix(x,eConn)
#=  Creates the 2D finite element mass matrix (Natural boundary conditions)

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  In the future, this should have an optional argument (ide) that allows
  for various boundary condition types.
=#

  nNodes    = size(x    ,1);
  nElements = size(eConn,1);

  #  Set the index into equation numbers
  ide = zeros(Int64,nNodes);

#  #  If there are no Dirichlet boundary conditions, and no ide array...
#  nUnk = 0;
#  iU = zeros(Int64,nNodes);
#  for i=1:nNodes
#    if ide[i]==0
#      nUnk = nUnk+1;
#      iU[nUnk] = i;
#      ide[i] = nUnk;
#    end
#  end

  ide = 1:nNodes;

  rule  = 3;
  r,s,w = twodQuadratureRule(rule);
  o     = ones(rule);

  nElDOF = size(eConn,2);
  nEntries = nElements*nElDOF^2;
  II  = zeros(Int64,nEntries);
  JJ  = zeros(Int64,nEntries);
  MM  = zeros(Float64,nEntries);

  xg  = zeros(Float64,rule,2);
  wg  = zeros(Float64,rule);
  phi = zeros(Float64,rule,nElDOF);
  p_x = zeros(Float64,rule,nElDOF);
  p_y = zeros(Float64,rule,nElDOF);

  index = 0;
  for k=1:nElements

    nLocal = eConn[k,:][:];
    xLocal = x[nLocal,:];

    xg,wg,phi,p_x,p_y = twodShape( xLocal, r, s, w )

    MLoc = twodBilinear(  o, phi, phi, wg )

    lDOF = ide[nLocal];
    for nt = 1:nElDOF
      for nu = 1:nElDOF
        index = index + 1;
        II[index] = lDOF[nu];
        JJ[index] = lDOF[nt];
        MM[index] = MLoc[nt,nu];
      end
    end
  end

  M = sparse(II,JJ,MM);

  return M
end


text:
using DrakeVisualizer
using GeometryTypes
using CoordinateTransformations
using Meshing
using Test
using Base.Iterators: product
using LinearAlgebra

include("lazytree.jl")
include("comms_t.jl")
include("visualizer.jl")
include("polyhedra.jl")
include("notebook.jl")


text:
"""
A file with helper functions for spatial simulations
"""

"""
stores info for a spatial jump
"""
struct SpatialJump{J}
    "source location"
    src::J

    "index of jump as a hop or reaction"
    jidx::Int

    "destination location, equal to src for within-site reactions"
    dst::J
end

function Base.show(io::IO, ::MIME"text/plain", jump::SpatialJump)
    println(io, "SpatialJump with source $(jump.src), destination $(jump.dst) and index $(jump.jidx).")
end

######################## helper routines for all spatial SSAs ########################
"""
    sample_jump_direct(p, site)

sample jump at site with direct method
"""
function sample_jump_direct(p, site)
    if rand(p.rng)*(total_site_rate(p.rx_rates, p.hop_rates, site)) < total_site_rx_rate(p.rx_rates, site)
        rx = sample_rx_at_site(p.rx_rates, site, p.rng)
        return SpatialJump(site, rx+p.numspecies, site)
    else
        species_to_diffuse, target_site = sample_hop_at_site(p.hop_rates, site, p.rng, p.spatial_system)
        return SpatialJump(site, species_to_diffuse, target_site)
    end
end

total_site_rate(rx_rates::RxRates, hop_rates::AbstractHopRates, site) = total_site_hop_rate(hop_rates, site) + total_site_rx_rate(rx_rates, site)

function update_rates_after_reaction!(p, u, site, reaction_id)
    update_rx_rates!(p.rx_rates, p.dep_gr[reaction_id], u, site)
    update_hop_rates!(p.hop_rates, p.jumptovars_map[reaction_id], u, site, p.spatial_system)
end

function update_rates_after_hop!(p, u, source_site, target_site, species)
    update_rx_rates!(p.rx_rates, p.vartojumps_map[species], u, source_site)
    update_hop_rate!(p.hop_rates, species, u, source_site, p.spatial_system)

    update_rx_rates!(p.rx_rates, p.vartojumps_map[species], u, target_site)
    update_hop_rate!(p.hop_rates, species, u, target_site, p.spatial_system)
end

"""
    update_state!(p, integrator)

updates state based on p.next_jump
"""
function update_state!(p, integrator)
    jump = p.next_jump
    if is_hop(p, jump)
        execute_hop!(integrator, jump.src, jump.dst, jump.jidx)
    else
        rx_index = reaction_id_from_jump(p,jump)
        @inbounds executerx!((@view integrator.u[:,jump.src]), rx_index, p.rx_rates.ma_jumps)
    end
    # save jump that was just exectued
    p.prev_jump = jump
    nothing
end

"""
    is_hop(p, jump)

true if jump is a hop
"""
function is_hop(p, jump)
    jump.jidx <= p.numspecies
end

"""
    execute_hop!(integrator, jump)

documentation
"""
function execute_hop!(integrator, source_site, target_site, species)
    @inbounds integrator.u[species,source_site] -= 1
    @inbounds integrator.u[species,target_site] += 1
end

"""
    reaction_id_from_jump(p,jump)

return reaction id by subtracting the number of hops
"""
function reaction_id_from_jump(p,jump)
    jump.jidx - p.numspecies
end


text:
@testset "Graph operations" begin
    g = SimpleGraph(5)
    @test_throws DomainError add_edges_distance!(g, 1.5; distance=5)
    @test_throws DomainError add_edges_distance!(g, 1; distance=-5)

    nb_edges_to_add = add_edges_distance!(g, 2; distance=2)
    @test ne(g) == nb_edges_to_add

    g = path_graph(5)
    nb_edges_to_add = add_edges_distance!(g, 2; distance=2)
    @test ne(g) - 4 == nb_edges_to_add
    added_edges = setdiff(edges(g), edges(path_graph(5)))
    for e in added_edges
        @test abs(e.src - e.dst) == 2
    end
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks