text:
using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, LinearAlgebra

@parameters x, t
@variables u(..)

Dx(d) = Differential(x)^d

Dt = Differential(t)

t_min= 0.
t_max = 2.0
x_min = 0.
x_max = 20.0

dx = 1.0

domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]

@testset "Test 01: Cartesian derivative discretization" begin
    weights = []
    push!(weights, ([-0.5,0,0.5], [1.,-2.,1.], [-1/2,1.,0.,-1.,1/2]))
    push!(weights, ([1/12, -2/3,0,2/3,-1/12], [-1/12,4/3,-5/2,4/3,-1/12], [1/8,-1.,13/8,0.,-13/8,1.,-1/8]))
    for d in 1:3
        for (i,a) in enumerate([2])
            pde = Dt(u(t,x)) ~ +Dx(d)(u(t,x))
            bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)]

            @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])

            # Test centered order
            disc = MOLFiniteDifference([x=>dx], t; approx_order=a)

            depvar_ops = map(x->operation(x.val),pdesys.depvars)

            depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
            depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
            depvars = collect(depvars_lhs ∪ depvars_rhs)
            # Read the independent variables,
            # ignore if the only argument is [t]
            indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
            x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))

            s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)

            derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)

            #@show pde.rhs, operation(pde.rhs), arguments(pde.rhs)
            for II in s.Igrid[s.ū[1]][2:end-1]
                #II = s.Igrid[end-1]
                I1 = MethodOfLines.unitindices(1)[1]

                rules = MethodOfLines.generate_finite_difference_rules(II, s, depvars, pde, derivweights)

                disc_pde=substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
                #@show disc_pde
                #@test disc_pde == dot(weights[i][d], s.discvars[depvars[1]][[II + j*I1 for j in MethodOfLines.half_range(length(weights[i][d]))]])
                ufunc(u, Itap, x) = s.discvars[u][Itap]
                if isodd(d)
                    @test isequal(substitute(pde.rhs,rules), MethodOfLines.upwind_difference(d, II, s, derivweights, (1, x), depvars[1], ufunc, true))
                else
                    @test isequal(substitute(pde.rhs,rules), MethodOfLines.central_difference(derivweights.map[Differential(x)^d], II, s, (1, x), depvars[1], ufunc))
                end
            end
        end
    end
end

@testset "Test 01: Nonlinear Diffusion discretization" begin

    pde = Dt(u(t,x)) ~ Dx(1)(u(t,x))
    bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])

    # Test centered order

    depvar_ops = map(x->operation(x.val),pdesys.depvars)

    depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
    depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
    depvars = collect(depvars_lhs ∪ depvars_rhs)
    # Read the independent variables,
    # ignore if the only argument is [t]
    indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
    x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))

    for order in [2]
        disc = MOLFiniteDifference([x=>dx], t; approx_order=order)
        s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)

        derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)

        ufunc(u, I, x) = s.discvars[u][I]
        #TODO Test Interpolation of params
        # Test simple case
        for II in s.Igrid[s.ū[1]][2:end-1]
            expr = MethodOfLines.cartesian_nonlinear_laplacian((1~1).lhs, II, derivweights, s, depvars, x, u(t,x))
            expr2 = MethodOfLines.central_difference(derivweights.map[Differential(x)^2], II, s, (1,x), u(t,x), ufunc)
            @test isequal(expr, expr2)
        end
    end
end

@testset "Test 02: Spherical Diffusion discretization" begin

    pde  = Dt(u(t,x)) ~ 1/x^2 * Dx(1)(x^2 * Dx(1)(u(t,x)))

    bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])

    # Test centered order
    disc = MOLFiniteDifference([x=>dx], t; approx_order=2)

    depvar_ops = map(x->operation(x.val),pdesys.depvars)

    depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
    depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
    depvars = collect(depvars_lhs ∪ depvars_rhs)
    # Read the independent variables,
    # ignore if the only argument is [t]
    indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
    x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))

    s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)

    derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)

    for II in s.Igrid[s.ū[1]][2:end-1]
        #TODO Test Interpolation of params
        expr = MethodOfLines.spherical_diffusion((1~1).lhs, II, derivweights, s, depvars, x, u(t,x))
        #@show II, expr
    end

end


text:
function handle_macro(@nospecialize(x), state) end
function handle_macro(x::EXPR, state)
    typof(x) !== MacroCall && return
    if typof(x.args[1]) === MacroName
        state(x.args[1])
        if _points_to_Base_macro(x.args[1], "deprecate", state) && length(x.args) == 3
            if bindingof(x.args[2]) !== nothing
                return
            elseif CSTParser.is_func_call(x.args[2])
                # add deprecated method
                # add deprecated function binding and args in new scope
                mark_binding!(x.args[2], x)
                mark_sig_args!(x.args[2])
                s0 = state.scope # store previous scope
                state.scope = Scope(s0, x, Dict(), nothing, false)
                setscope!(x, state.scope) # tag new scope to generating expression
                state(x.args[2])
                state(x.args[3])
                state.scope = s0
            elseif isidentifier(x.args[2])
                mark_binding!(x.args[2], x)
            end
        elseif _points_to_Base_macro(x.args[1], "enum", state)
            for i = 2:length(x.args)
                if !(typof(x.args[i]) === PUNCTUATION)
                    if bindingof(x.args[i]) !== nothing
                        break
                    end
                    if i == 3 && typof(x.args[i]) === CSTParser.Begin
                        for j in 1:length(x.args[i].args[2].args)
                            mark_binding!(x.args[i].args[2].args[j], x)
                        end
                        break
                    end
                    mark_binding!(x.args[i], x)
                end
            end
        elseif _points_to_Base_macro(x.args[1], "goto", state)
            if length(x.args) == 2 && typof(x.args[2]) === CSTParser.IDENTIFIER
                setref!(x.args[2], Binding(noname, nothing, nothing, EXPR[], nothing, nothing))
            end
        elseif _points_to_Base_macro(x.args[1], "label", state)
            if length(x.args) == 2 && typof(x.args[2]) === CSTParser.IDENTIFIER
                mark_binding!(x.args[2])
            end
        # elseif _points_to_Base_macro(x.args[1], "nospecialize", state)
        elseif length(x.args[1].args) == 2 && isidentifier(x.args[1].args[2]) && valof(x.args[1].args[2]) == "nospecialize"
            for i = 2:length(x.args)
                if !(typof(x.args[i]) === PUNCTUATION)
                    if bindingof(x.args[i]) !== nothing
                        break
                    end
                    mark_binding!(x.args[i], x)
                end
            end
        elseif _points_to_arbitrary_macro(x.args[1], "Turing", "model", state) && length(x.args) == 2 &&
            _binary_assert(x.args[2], CSTParser.Tokens.EQ) &&
            _expr_assert(x.args[2].args[3], CSTParser.Begin, 3) && typof(x.args[2].args[3].args[2]) === CSTParser.Block && x.args[2].args[3].args[2].args isa Vector{EXPR}
            for i = 1:length(x[2][3][2].args)
                ex = x[2][3][2][i]
                if typof(ex) == CSTParser.BinaryOpCall && kindof(ex[2]) === CSTParser.Tokens.APPROX
                    mark_binding!(ex)
                end
            end
        elseif _points_to_arbitrary_macro(x.args[1], "JuMP", "variable", state)
            if length(x.args) < 3
                return
            elseif length(x.args) >= 5 && typof(x.args[2]) === PUNCTUATION
                _mark_JuMP_binding(x.args[5])
            else
                _mark_JuMP_binding(x.args[3])
            end
        elseif (_points_to_arbitrary_macro(x.args[1], "JuMP", "expression", state) ||
            _points_to_arbitrary_macro(x.args[1], "JuMP", "NLexpression", state) ||
            _points_to_arbitrary_macro(x.args[1], "JuMP", "constraint", state) || _points_to_arbitrary_macro(x.args[1], "JuMP", "NLconstraint", state)) && length(x.args) > 1
            if typof(x.args[2]) === PUNCTUATION
                if length(x.args) == 8
                    _mark_JuMP_binding(x.args[5])
                end
            else
                if length(x.args) == 4
                    _mark_JuMP_binding(x.args[3])
                end
            end
        end
    elseif typof(x.args[1]) == CSTParser.GlobalRefDoc && length(x.args) == 3 && CSTParser.isliteral(x.args[2]) && kindof(x.args[2]) === CSTParser.Tokens.TRIPLE_STRING && isidentifier(x.args[3])
        mark_binding!(x.args[3])
        setref!(x.args[3], bindingof(x.args[3]))
    end
end

function _rem_ref(x::EXPR)
    if typof(x) === CSTParser.Ref && x.args isa Vector{EXPR} && length(x.args) > 0
        return x.args[1]
    end
    return x
end

function _mark_JuMP_binding(arg)
    if CSTParser.isidentifier(arg) || typof(arg) === CSTParser.Ref
        mark_binding!(_rem_ref(arg))
    elseif _binary_assert(arg, CSTParser.Tokens.EQEQ) || _binary_assert(arg, CSTParser.Tokens.LESS_EQ)  || _binary_assert(arg, CSTParser.Tokens.GREATER_EQ)
        if CSTParser.isidentifier(arg.args[1]) || typof(arg.args[1]) === CSTParser.Ref
            mark_binding!(_rem_ref(arg.args[1]))
        else
            mark_binding!(_rem_ref(arg.args[3]))
        end
    elseif typof(arg) === CSTParser.Comparison && arg.args isa Vector{EXPR} && length(arg.args) == 5
        mark_binding!(_rem_ref(arg.args[3]))
    end
end

function _points_to_Base_macro(x::EXPR, name, state)
    length(x.args) == 2 && isidentifier(x.args[2]) && valof(x.args[2]) == name && refof(x.args[2]) == getsymbolserver(state.server)["Base"].vals[string("@", name)]
end

function _points_to_arbitrary_macro(x::EXPR, module_name, name, state)
    length(x.args) == 2 && isidentifier(x.args[2]) && valof(x.args[2]) == name && haskey(getsymbolserver(state.server), module_name) && haskey(getsymbolserver(state.server)[module_name].vals, string("@", name)) && refof(x.args[2]) == getsymbolserver(state.server)[module_name].vals[string("@", name)]
end

ext: jl


text:
using DynamicalSystems
using OrdinaryDiffEq
using GLMakie
using DataStructures: CircularBuffer

# Function defining ODEs for model
function model!(du, u, p, t)
	ω0, ωF, ωM, ϕF, ϕM, μ, ν, αF, αM, ψF, ψM = p
	du[1] = ω0[]*π
	du[2] = ωF[]*π + μ[]*sin(u[3]-u[2]-ϕF[]*π) + αF[]*sin(u[1]-u[2]-ψF[]*π)
	du[3] = ωM[]*π + ν[]*sin(u[2]-u[3]-ϕM[]*π) + αM[]*sin(u[1]-u[3]-ψM[]*π)
end

# Function to convert phase states to x,y coordinates
function xycoords(state)
    x1 = sin(state[1])
    y1 = cos(state[1])
    x2 = sin(state[2])
    y2 = cos(state[2])
	x3 = sin(state[3])
	y3 = cos(state[3])
    return x1,x2,x3,y1,y2,y3
end

# Function to iterate system state
function progress_for_one_step!(integ)
    step!(integ)
    u = integ.u
    return xycoords(u)
end

# Function to update figure based on system iteration
function animstep!(integ, dots, hands)
    x1,x2,x3,y1,y2,y3 = progress_for_one_step!(integ)
    dots[] = [Point2f0(x1, y1), Point2f0(x2, y2), Point2f0(x3, y3)]
	hands[] = [Point2f0(x1, y1), Point2f0(0,0), Point2f0(x2, y2), Point2f0(0, 0), Point2f0(x3,y3)]
end

# Set up figure canvas
fig = Figure(); display(fig)
ax = Axis(fig[1,1])
ax.title = "Circadian Rhythms"
ax.aspect = DataAspect() # ??
lim = 1.1 # Plot canvas limit
xlims!(ax, -lim, lim)
ylims!(ax, -lim, lim)
lines!(ax,[Point2f0(0.0,lim),Point2f0(0.0,-lim)]; linestyle=:dash, alpha=0.5, color=:black)
lines!(ax,[Point2f0(lim,0.0),Point2f0(-lim,0.0)]; linestyle=:dash, alpha=0.5, color=:black)
hidedecorations!(ax)

# Initial conditions
u0   = rand(3).*(2π) # Initial phases (Sunlight, fibroblast, macrophage)
# Plot initial conditions as a vector of observable points.
x1,x2,x3,y1,y2,y3 = xycoords(u0)
# 3 points: one for each cell type
dots = Observable([Point2f0(x1, y1), Point2f0(x2, y2), Point2f0(x3, y3)])
# 3 lines corresponding to a clock hand for each cell type
hands= Observable([Point2f0(x1, y1), Point2f0(0,0), Point2f0(x2, y2), Point2f0(0, 0), Point2f0(x3,y3)])
# Plot objects
lines!(ax, hands; linewidth = 4, color = :black)
scatter!(ax,dots;marker=[:star8,:circle,:circle],color=[:red,:green,:blue],markersize=48)

# Set up parameter sliders
# ω0  Sunlight phase rate of change /day
# ωF  Fibroblast phase intrinsic rate of change /day
# ωM = 2π*24/23.5              # Macrophage phase intrinsic rate of change /day
# ϕF = 0.0            # Fibroblast phase offset in coupling to macrophage
# ϕM = 0.0 		     # Macrophage phase offset in coupling to fibroblast
# μ  = 10.0             # Amplitude of fibroblast phase coupling to macrophage
# ν  = 1.0             # Amplitude of macrophage phase coupling to fibroblast
# αF = 0.0             # Amplitude of fibroblast phase coupling to sunlight
# αM = 1.0             # Amplitude of macrophage phase coupling to sunlight
# ψF = 0.0             # Fibroblast phase offset from sunlight
# ψM = 0.0             # Macrophage phase offset from sunlight
lsgrid = labelslidergrid!(
    fig,
    ["ω0", "ωF", "ωM", "ϕF", "ϕM", "μ", "ν", "αF", "αM", "ψF", "ψM"],
    [0:0.1:2, 0:0.1:2, 0:0.1:2, 0:0.1:2, 0:0.1:2, 0:0.1:5.0, 0:0.1:5.0, 0:0.1:5.0, 0:0.1:5.0, 0:0.1:2, 0:0.1:2];
    formats = [x -> "$(round(x, digits = 1))$s" for s in ["π", "π", "π", "π", "π", "", "", "", "", "π", "π"]],
    width = 350,
    tellheight = false
)
fig[1, 2] = lsgrid.layout
# Set default slider values
defaults = [2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
set_close_to!.(lsgrid.sliders, defaults)

# Pull parameters from slider positions [ω0, ωF, ωM, ϕF, ϕM, μ, ν, αF, αM, ψF, ψM]
sliderobservables = [s.value for s in lsgrid.sliders]

# Set up differential equation integrator
prob = ODEProblem(model!,u0,(0.0,10000.0),sliderobservables)
dp = ContinuousDynamicalSystem(prob)
# Solve diffeq with constant step for smoother curves
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.005)
# Set up integrator for each iteration
integ = integrator(dp, u0; diffeq...)

run = Button(fig[2,1]; label = "Start/Stop", tellwidth = false)
isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
    @async while isrunning[]
        isopen(fig.scene) || break # ensures computations stop if closed window
        animstep!(integ, dots, hands)
        sleep(0.02) # or `yield()` instead
    end
end

max_stars_repo_name: chris-revell/ClockOscillators

max_issues_repo_name: chris-revell/ClockOscillators

max_forks_repo_name: chris-revell/ClockOscillators


text:
"""
温度流的Bubble
"""

"""
获取所有的Bubble(pp, fs, ex)
qpp = k1+k2; qfs = k3-k2; qex = k1-k3
"""
function all_bubble_tf_mt(Γ4::Gamma4{T, P}, lval) where {T, P}
    #
    lamb = Γ4.λ_0 * exp(-lval)
    brlu_area = area(Γ4.model.brillouin)
    #因为现在所有的能带都必须有同一个lpats，所以只要每一个patch中有哪些tri
    tris_pat = Γ4.ltris_pat
    #
    bubbval_pp = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b2
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b1), k2(b2)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_pp)
        alpha, beta, b1, b2, i_n, i1, i2 = Tuple(idxs)
        k1, k2 = Γ4.patches[b1][i1], Γ4.patches[b2][i2]
        q_pp = Γ4.model.kadd(k1, k2)
        bubbres = pi_αβ_minus_tf(
            tris_pat[i_n], brlu_area, lamb,
            q_pp, Γ4.model.dispersion[alpha], Γ4.model.dispersion[beta]
        )
        bubbval_pp[alpha, beta, b1, b2, i_n, i1, i2] = bubbres
    end
    bubb_qpp = Bubble{:tf, :minus}(lval, bubbval_pp)
    #
    bubbval_fs = Array{Float64, 7}(
        undef,
        #alpha, beta, b2, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k2(b2), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_fs)
        alpha, beta, b2, b3, i_n, i2, i3 = Tuple(idxs)
        k2, k3 = Γ4.patches[b2][i2], Γ4.patches[b3][i3]
        q_fs = Γ4.model.kadd(k3, -k2)
        bubbres = pi_αβ_plus_tf(
            tris_pat[i_n], brlu_area,
            lamb, q_fs, Γ4.model.dispersion[alpha],
            Γ4.model.dispersion[beta]
        )
        bubbval_fs[alpha, beta, b2, b3, i_n, i2, i3] = bubbres
    end
    bubb_qfs = Bubble{:tf, :plus}(lval, bubbval_fs)
    #
    bubbval_ex = Array{Float64, 7}(
        undef,
        #alpha, beta, b1, b3
        Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum, Γ4.model.bandnum,
        #nidx(alpha), k1(b3), k3(b3)
        Γ4.patchnum, Γ4.patchnum, Γ4.patchnum
    )
    Threads.@threads for idxs in CartesianIndices(bubbval_ex)
        alpha, beta, b1, b3, i_n, i1, i3 = Tuple(idxs)
        k1, k3 = Γ4.patches[b1][i1], Γ4.patches[b3][i3]
        q_ex = Γ4.model.kadd(k1, -k3)
        bubbres = pi_αβ_plus_tf(
            tris_pat[i_n], brlu_area, lamb,
            q_ex, Γ4.model.dispersion[alpha],
            Γ4.model.dispersion[beta]
        )
        bubbval_ex[alpha, beta, b1, b3, i_n, i1, i3] = bubbres
    end
    bubb_qex = Bubble{:tf, :plus}(lval, bubbval_ex)
    return bubb_qpp, bubb_qfs, bubb_qex
end

"""
温度流的Bubble
"""
function pi_αβ_plus_tf(
    ltris::Vector{T},
    area::Float64, lamb::Float64,
    qval::Point2D,
    dispα::Function, dispβ::Function) where T <: Basics.AbstractTriangle
    """温度流的+
    这里的lamb就是T，ltris中的所有三角都应该要在同一个patch中,
    tarea是每个小三角形的面积(已经不用了)，
    dispa是和k相关的那个能带，dispb是k-q相关的
    """
    nega_q = -qval
    result = 0.
    for tri in ltris
        #这个小三角形的k值
        kval = tri.center
        #k-q
        kprim = kval + nega_q
        #kprim = kadd(kval, nega_q)
        #epsilon_k
        eps_k = dispα(kval.x, kval.y)
        #epsilon_{k-q}
        eps_kp = dispβ(kprim.x, kprim.y)
        #这个小区域的贡献
        if abs(eps_k - eps_kp) < 1.e-10
            #如果特别小，可以利用
            # lim (eps_k -> eps_kp) Pi^{+} =
            # 1/T (e^{eps/T} (-eps/T*e^{eps/T} +
            #eps/T + e^{eps/T} + 1)) / (e^{eps/T} + 1)^3
            bval = eps_kp / lamb
            #如果本身就很大，分母会比较大导致接近0
            if bval > 25
                #@info "数值不稳定"
                num = 0.
                den = 100.
            else
                expb = exp(bval)
                num = expb * (-bval * expb + bval + expb + 1)
                den = (1+expb)^3
            end
            d_val = num / den / lamb
        else#如果本身是正常的数值
            #分子左侧的数值
            if (eps_k / lamb) > 25
                #@info "数值不稳定"
                num_left = 0.
            else
                #exp^{epsilon_k / T}
                exp_k_t = exp(eps_k / lamb)
                num_left = eps_k / lamb * exp_k_t / ((1 + exp_k_t)^2)
            end
            #分子右侧的数值
            if (eps_kp / lamb) > 25
                #@info "数值不稳定"
                num_righ = 0.
            else
                #e^{epsilon_{k-q} / T}
                exp_kp_t = exp(eps_kp / lamb)
                num_righ = eps_kp / lamb * exp_kp_t / ((1 + exp_kp_t)^2)
            end
            d_val = (num_left - num_righ) / (eps_k - eps_kp)
        end# end if 小区域的贡献计算完
        result += d_val * tri.area
    end#对ltris的循环
    result = result / area
    return result
end

function pi_αβ_minus_tf(
    ltris::Vector{T},
    area::Float64, lamb::Float64,
    qval::Point2D,
    dispα::Function, dispβ::Function) where T <: Basics.AbstractTriangle
    result = 0.
    for tri in ltris
        #这个小三角形的k值
        kval = tri.center
        nega_k = -kval
        #-k+q
        kprim = nega_k + qval
        #kprim = kadd(nega_k, qval)
        #epsilon_k
        eps_k = dispα(kval.x, kval.y)
        #-epsilon_{-k+q}
        neps_kp = -dispβ(kprim.x, kprim.y)
        #这个时候，因为epsilon_{-k+q}前面已经有了负号，分母上还是负号
        #计算这个小区域的贡献
        if abs(eps_k - neps_kp) < 1.e-10
            #如果两个数值比较接近, Pi^{-}和Pi^{+}的公式完全一样，就是第二个能量要加个负号
            # lim (eps_k -> -eps_kp) Pi^{-} =
            # 1/T (e^{eps/T} (-eps/T*e^{eps/T} + eps/T + e^{eps/T} + 1))
            #/ (e^{eps/T} + 1)^3
            bval = eps_k / lamb
            if bval > 25
                #@info "数值不稳定"
                num = 0.
                den = 100.
            else
                expb = exp(bval)
                num = expb * (-bval * expb + bval + expb + 1)
                den = (1+expb)^3
            end
            d_val = num / den / lamb
        else
            if (eps_k / lamb) > 25
                #@info "数值不稳定"
                num_left = 0.
            else
                #e^{epsilon_k / T}
                exp_k_t = exp(eps_k / lamb)
                num_left = eps_k / lamb * exp_k_t / ((1 + exp_k_t)^2)
            end
            if (neps_kp / lamb) > 25
                #@info "数值不稳定"
                num_righ = 0.
            else
                #e^{-epsilon_{-k+q} / T}
                exp_nkp_t = exp(neps_kp / lamb)
                num_righ = neps_kp / lamb * exp_nkp_t /
                ((1 + exp_nkp_t)^2)
            end
            d_val = (num_left - num_righ) / (eps_k - neps_kp)
        end # end if 小区域的贡献计算完
        result += d_val * tri.area
    end#对ltris的循环
    result = result / area
    return result
end


text:
using Documenter, MethodOfLines

# Make sure that plots don't throw a bunch of warnings / errors!
ENV["GKSwstype"] = "100"

include("pages.jl")

makedocs(
    sitename="MethodOfLines.jl",
    authors="Chris Rackauckas, Alex Jones et al.",
    clean=true,
    doctest=false,
    strict=[
        :doctest,
        :linkcheck,
        :parse_error,
        :example_block,
        # Other available options are
        # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
    ],
    modules=[MethodOfLines], format=Documenter.HTML(analytics="UA-90474609-3",
        assets=["assets/favicon.ico"],
        canonical="https://methodoflines.sciml.ai/stable/"), pages=pages
)

deploydocs(
    repo="github.com/SciML/MethodOfLines.jl";
    push_preview=true
)


text:
include("../utils.jl")

@testset "Scatter2D" begin
    d = [1 1; 2 2]
    s = G.Scatter2D(d, false, 1)
    @test s isa G.Drawing2D
    @test s isa G.Drawing
    @test s isa G.Scatter2D{Matrix{T}} where T
    @test s.hasmissing == false
    @test s.nobj == 1
    @test length(s.linestyles) == 1
    @test length(s.markerstyles) == 1
    @test isempty(s.labels)

    dh = G.DrawingHandle(s)
    @test dh isa G.DrawingHandle{G.Scatter2D{T}} where T
end

@testset "Fill2D" begin
    d = [1 1; 2 2]
    s = G.Fill2D(d, 0.5)
    @test s isa G.Fill2D{Matrix{T}} where T
    @test s.xmin == 0.5
    @test !G.isdef(s.xmax)
    @test isempty(s.label)
end

@testset "Hist2D" begin
    d = [1 1; 2 2]
    s = G.Hist2D(d, false, 3, (0.5,0.7))
    @test s isa G.Hist2D{Matrix{T}} where T
    @test s.hasmissing == false
    @test s.nobs == 3
    @test s.range == (0.5, 0.7)

    @test s.horiz == false
    @test s.bins == nothing
    @test s.scaling == "none"
    @test s.label == ""
end

@testset "Bar2D" begin
    d = [1 1; 2 2]
    s = G.Bar2D(d, false, 2)
    @test s isa G.Bar2D{Matrix{T}} where T
    @test length(s.barstyles) == 2
    @test s.stacked == false
    @test s.horiz == false
    @test s.bwidth === nothing
    @test isempty(s.labels)
end

@testset "Boxplot" begin
    d = [0.1 0.2; 0.3 0.4]
    b = G.Boxplot(d, 2)
    @test b isa G.Boxplot
    @test length(b.boxstyles) == b.nobj == 2
    @test b.horiz == false
end

@testset "Heatmap" begin
    d = [1 2; 1 2]
    h = G.Heatmap(d, 0.1, 0.2)
    @test h isa G.Heatmap
    @test h.zmin == 0.1
    @test h.zmax == 0.2
    @test length(h.cmap) == 10
    @test h.cmiss == c"white"
    @test h.transpose == false
end


text:
using Flux.Tracker
using Flux

bandits = [0.2,0,-0.2,-5]

num_bandits = length(bandits)

function pullBandit(bandit)
    result = rand()
    if result > bandit
        return 1
    else
        return -1
    end
end

W = param(ones(num_bandits))

function reward(a)
    pullBandit(bandits[a])
end

λ = 0.001
ϵ = 0.1
loss() = -(log(W[chosen_action])*reward(chosen_action))
chosen_action = rand(1:4)

for i in 1:100000
    chosen_action = rand() < ϵ ? rand(1:num_bandits) : findmax(W.data)[2]
    back!(loss())
    W.data[chosen_action] -= λ*(W.grad[chosen_action])
    W.grad .= 0.0
end

W

max_stars_repo_name: fisimon/RL-for-FIFA-world-cup

max_issues_repo_name: fisimon/RL-for-FIFA-world-cup

max_forks_repo_name: fisimon/RL-for-FIFA-world-cup


text:
@testset "Gradient" begin
    @testset "Gradient in physical space" begin
        p = Parameter(2.0, :p)
        x = RandomVariable(Uniform(-2.0, 2.0), :x)
        y = RandomVariable(Uniform(-2.0, 2.0), :y)

        model = Model(df -> df.p .* df.x .^ 2 - df.y .^ 2, :f)

        g = gradient([model], [p, x, y], DataFrame(; p=2.0, x=2.0, y=1.0), :f)

        @test isapprox(g.x, 8.0; rtol=0.001)
        @test isapprox(g.y, -2.0; rtol=0.001)
    end
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
function _get_organism_in_front(model::GebModel, organism::GebOrganism)
    _get_organism_in_direction(model, organism.coordinates, organism.direction)
end
function _get_organism_behind(model::GebModel, organism::GebOrganism)
    _get_organism_in_direction(model, organism.coordinates, 360 - organism.direction)
end

function _get_organism_in_direction(model::GebModel, coordinates, direction)
    x_target, y_target = _get_cell_in_direction(model, coordinates, direction)
    target = model.grid[x_target, y_target]
    return target
end

function _get_cell_in_direction(model, coordinates, direction)
    x_organism, y_organism = _get_grid_coordinates(coordinates)

    if 45 <= direction < 135
        x_target, y_target = x_organism, y_organism + 1
    elseif 135 <= direction < 225
        x_target, y_target = x_organism - 1, y_organism
    elseif 225 <= direction < 315
        x_target, y_target = x_organism, y_organism - 1
    elseif 315 <= direction || direction < 45
        x_target, y_target = x_organism + 1, y_organism
    end

    x_target, y_target = mod_1(x_target, model.size), mod_1(y_target, model.size)

    return x_target, y_target
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text: distance(x, y) = sqrt(sum((x - y).^2))


text:
include("linesearch/ls_approach.jl")
include("trustregion/tr_approach.jl")

function trace_show(show_trace, f_curr, f_next, x_next, x_curr, α)
    if show_trace
        println("Objective value (curr): ", f_next)
        println("Objective value (prev): ", f_curr)
        println("Step size: ", α)
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Bitmap/Histogram

using Color, Images, FixedPointNumbers

ima = imread("bitmap_histogram_in.jpg")
imb = convert(Image{Gray{Ufixed8}}, ima)

# calculate histogram
a = map(x->x.val.i, imb.data)
(nothing, h) = hist(reshape(a, length(a)), -1:typemax(Uint8))

g = float(imb.data)
b = g .> median(g)
fill!(imb, Gray(0.0))
imb[b] = Gray(1.0)

imwrite(imb, "bitmap_histogram_out.png")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
export VRepSolver, VRepPolyhedraModel

"""
    VRepSolver

Linear Programming solver using the V-representation of the feasible set to find the optimal solution.
"""
struct VRepSolver <: MPB.AbstractMathProgSolver
end

mutable struct VRepPolyhedraModel <: AbstractPolyhedraModel
    vrep::Union{VRep, Nothing}
    obj::Union{Vector, Nothing}
    sense::Symbol

    objval
    solution::Union{AbstractVector, Nothing}
    status::Symbol

    function VRepPolyhedraModel()
        new(nothing, nothing, :Feas, nothing, nothing, :Undecided)
    end
end

PolyhedraModel(::VRepSolver) = VRepPolyhedraModel()
MPBSI.LinearQuadraticModel(s::VRepSolver) = PolyhedraToLPQPBridge(PolyhedraModel(s))

function MPBSI.loadproblem!(lpm::VRepPolyhedraModel, vrep::VRep, obj, sense)
    if !(sense in [:Max, :Min])
        error("sense should be :Max or :Min")
    end
    lpm.vrep = vrep
    if isapproxzero(obj)
        lpm.sense = :Feas
    else
        lpm.obj = copy(obj)
        lpm.sense = sense
    end
end

function MPBSI.optimize!(lpm::VRepPolyhedraModel)
    if lpm.vrep === nothing
        error("No problem loaded")
    end
    prob = lpm.vrep
    N = fulldim(prob)
    T = MultivariatePolynomials.coefficienttype(prob)
    lpm.status = :Undecided
    if !haspoints(prob) && !haslines(prob) && !hasrays(prob)
        lpm.status = :Infeasible
    else
        if lpm.sense in [:Max, :Min]
            better(a, b) = (lpm.sense == :Max ? a > b : a < b)
            _better(a, b) = (lpm.sense == :Max ? _gt(a, b) : _lt(a, b))
            bestobjval = zero(T)
            for r in allrays(prob)
                objval = lpm.obj ⋅ r
                if _better(objval, bestobjval)
                    bestobjval = objval
                    lpm.solution = coord(r)
                end
            end
            if _better(bestobjval, zero(T))
                lpm.status = :Unbounded
                lpm.objval = lpm.sense == :Max ? typemax(T) : typemin(T)
            end
            if lpm.status != :Unbounded
                for p in points(prob)
                    objval = lpm.obj ⋅ p
                    if lpm.status == :Undecided || better(objval, lpm.objval)
                        lpm.status = :Optimal
                        lpm.objval = objval
                        lpm.solution = p
                    end
                end
            end
        else
            lpm.status = :Optimal
            lpm.objval = zero(T)
            lpm.solution = spzeros(T,N)
        end
    end
end

function MPB.status(lpm::VRepPolyhedraModel)
    lpm.status
end
function MPB.getobjval(lpm::VRepPolyhedraModel)
    lpm.objval
end
function MPB.getsolution(lpm::VRepPolyhedraModel)
    copy(lpm.solution)
end
function MPB.getunboundedray(lpm::VRepPolyhedraModel)
    copy(lpm.solution)
end


text:
using ReachabilityAnalysis: reach_homog_QINT

@testset "QINT algorithm" begin
    # x' = x - x^2, with solution x(t) = x0*e^t / (1 - (1 - e^t)x0)
    sol = reach_homog_QINT(a=-1.0, b=1.0, c=0.0, X0=Interval(0.1, 0.3),
                           T=10.0, Δ=0.1, δ=0.005, θ=0.01, maxiter=10_000)
    @test ρ([1.0], sol) < 1.02
end


text:
using Pkg;
Pkg.add("Plots")
Pkg.add("LinearAlgebra")
Pkg.add("Statistics")
Pkg.add("BenchmarkTools")
Pkg.add("Interpolations")
Pkg.add("Parameters")
Pkg.add("QuantEcon")
Pkg.add("Roots")
Pkg.add("Optim")
Pkg.add("Random")
using LinearAlgebra, Statistics
using BenchmarkTools, Interpolations, Parameters, Plots, QuantEcon, Roots
using Optim, Random
gr(fmt = :png);

function K!(Kg, g, grid, β, ∂u∂c, f, f′, shocks)
# This function requires the container of the output value as argument Kg

    # Construct linear interpolation object
    g_func = LinearInterpolation(grid, g, extrapolation_bc = Line())

    # solve for updated consumption value
    for (i, y) in enumerate(grid)
        function h(c)
            vals = ∂u∂c.(g_func.(f(y - c) * shocks)) .* f′(y - c) .* shocks
            return ∂u∂c(c) - β * mean(vals)
        end
        Kg[i] = find_zero(h, (1e-10, y - 1e-10))
    end
    return Kg
end

# The following function does NOT require the container of the output value as argument
K(g, grid, β, ∂u∂c, f, f′, shocks) =
    K!(similar(g), g, grid, β, ∂u∂c, f, f′, shocks)

function T(w, grid, β, u, f, shocks, Tw = similar(w);
		compute_policy = false)

	# apply linear interpolation to w
	w_func = LinearInterpolation(grid, w, extrapolation_bc = Line())

	if compute_policy
		σ = similar(w)
	end

	# set Tw[i]  = max_c { u(c) + β E w(f(y - c) z) }
	for (i, y) in enumerate(grid)
		objective(c) = u(c) + β * mean(w_func.(f(y - c) .* shocks))
		res = maximize(objective, 1e-10, y)

		if compute_policy
			σ[i] = Optim.maximizer(res)
		end
		Tw[i] = Optim.maximum(res)
	end

	if compute_policy
		return Tw, σ
	else
		return Tw
	end
end

isoelastic(c, γ) = isone(γ) ? log(c) : (c^(1 - γ) - 1) / (1 - γ)
Model  = @with_kw (α = 0.65,				#Productivity parameter
		   β = 0.95,				#Discount factor
		   γ = 1.0,				#Risk aversion
		   μ = 0.0,				#First parameter in lognorm(μ, σ)
		   s = 0.1,				#second parameter in lognorm
		   grid = range(1e-6, 4, length = 200),	#Grid
		   grid_min = 1e-6,			#Smallest grid point
		   grid_max = 4.0,			#Largest grid point
		   grid_size = 200,			#Number of grid points
		   u = (c, γ = γ) -> isoelastic(c, γ),	#utility function
		   ∂u∂c = c -> c^(-γ),			#u'
		   f = k -> k^α,			#production function
		   f' = k -> α * k^(α - 1)		#f'
		   )

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
# mathematical constants related to statistics

@irrational loghalf -0.6931471805599453094 log(big(0.5))
@irrational logtwo 0.6931471805599453094 log(big(2.))
@irrational logπ   1.1447298858494001741 log(big(π))
@irrational log2π  1.8378770664093454836 log(big(2.)*π)
@irrational log4π  2.5310242469692907930 log(big(4.)*π)


text:
"""
Subroutine dtrqsol

This subroutine computes the largest (non-negative) solution
of the quadratic trust region equation

  ||x + sigma*p|| = delta.

The code is only guaranteed to produce a non-negative solution
if ||x|| <= delta, and p != 0. If the trust region equation has
no solution, sigma = 0.

MINPACK-2 Project. March 1999.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dtrqsol(n,x,p,delta)
    zero = 0.0
    sigma = zero

    ptx = ddot(n,p,1,x,1)
    ptp = ddot(n,p,1,p,1)
    xtx = ddot(n,x,1,x,1)
    dsq = delta^2

    # Guard against abnormal cases.
    rad = ptx^2 + ptp*(dsq - xtx)
    rad = sqrt(max(rad,zero))

    if ptx > zero
        sigma = (dsq - xtx)/(ptx + rad)
    elseif rad > zero
        sigma = (rad - ptx)/ptp
    else
        sigma = zero
    end

    return sigma
end


text:
using Photometry
using Test

@testset "Photometry.jl" begin
    include("overlap.jl")
    include("circular.jl")
    include("photometry.jl")
end


text:
export chbv, chbv!

const α0 = 0.183216998528140087E-11

#=
Negative of the coefficients and poles of the partial fraction expansion of the
uniform rational Chebyshev approximation of type (14, 14) to the exponential
function over the negative real axis.
=#
const α = [-0.557503973136501826E+02 + im*0.204295038779771857E+03,
            0.938666838877006739E+02 - im*0.912874896775456363E+02,
           -0.469965415550370835E+02 + im*0.116167609985818103E+02,
            0.961424200626061065E+01 + im*0.264195613880262669E+01,
           -0.752722063978321642E+00 - im*0.670367365566377770E+00,
            0.188781253158648576E-01 + im*0.343696176445802414E-01,
           -0.143086431411801849E-03 - im*0.287221133228814096E-03]

const θ = [0.562314417475317895E+01 - im*0.119406921611247440E+01,
           0.508934679728216110E+01 - im*0.358882439228376881E+01,
           0.399337136365302569E+01 - im*0.600483209099604664E+01,
           0.226978543095856366E+01 - im*0.846173881758693369E+01,
          -0.208756929753827868E+00 - im*0.109912615662209418E+02,
          -0.370327340957595652E+01 - im*0.136563731924991884E+02,
          -0.889777151877331107E+01 - im*0.166309842834712071E+02]

const αconj = [-0.557503973136501826E+02 - im*0.204295038779771857E+03,
                0.938666838877006739E+02 + im*0.912874896775456363E+02,
               -0.469965415550370835E+02 - im*0.116167609985818103E+02,
                0.961424200626061065E+01 - im*0.264195613880262669E+01,
               -0.752722063978321642E+00 + im*0.670367365566377770E+00,
                0.188781253158648576E-01 - im*0.343696176445802414E-01,
               -0.143086431411801849E-03 + im*0.287221133228814096E-03]

const θconj = [0.562314417475317895E+01 + im*0.119406921611247440E+01,
               0.508934679728216110E+01 + im*0.358882439228376881E+01,
               0.399337136365302569E+01 + im*0.600483209099604664E+01,
               0.226978543095856366E+01 + im*0.846173881758693369E+01,
              -0.208756929753827868E+00 + im*0.109912615662209418E+02,
              -0.370327340957595652E+01 + im*0.136563731924991884E+02,
              -0.889777151877331107E+01 + im*0.166309842834712071E+02]

"""
    chbv(A, vec)

Calculate matrix exponential acting on some vector using the Chebyshev method.

# Input

- `A`   -- matrix which can be dense or sparse
- `vec` -- vector on which the matrix exponential of `A` is applied

# Notes

This Julia implementation is based on Expokit's CHBV Matlab code by
Roger B. Sidje, see below.

---

    y = chbv( H, x )
    CHBV computes the direct action of the matrix exponential on
    a vector: y = exp(H) * x. It uses the partial fraction expansion of
    the uniform rational Chebyshev approximation of type (14,14).
    About 14-digit accuracy is expected if the matrix H is symmetric
    negative definite. The algorithm may behave poorly otherwise.
    See also PADM, EXPOKIT.

    Roger B. Sidje (rbs@maths.uq.edu.au)
    EXPOKIT: Software Package for Computing Matrix Exponentials.
    ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
"""
function chbv(A, vec::Vector{T}) where {T}
    result = convert(Vector{promote_type(eltype(A), T)}, copy(vec))
    return chbv!(result, A, vec)
end

chbv!(A, vec::Vector{T}) where {T} = chbv!(vec, A, copy(vec))

function chbv!(w::Vector{T}, A, vec::Vector{T}) where {T<:Real}
    p = min(length(θ), length(α))
    rmul!(copyto!(w, vec), α0)
    @inbounds for i = 1:p
        w .+= real((A - θ[i]*I) \ (α[i] * vec))
    end
    return w
end

function chbv!(w::Vector{T}, A, vec::Vector{T}) where {T<:Complex}
    p = min(length(θ), length(α))
    rmul!(copyto!(w, vec), α0)
    t = [θ; θconj]
    a = 0.5 * [α; αconj]
    @inbounds for i = 1:2*p
        w .+= (A - t[i]*I) \ (a[i] * vec)
    end
    return w
end # chbv!


text:
mutable struct RowTable <: AbstractRowTable
    rows::Vector
    colindex::CIndex
    function RowTable(rows::Vector{Any}, colindex::CIndex)
        new(rows, colindex)
    end
end

# For efficiency, we either need to use these methods, or else
# make a typed RowTable
# RowTable(rows::Vector, colindex::CIndex) = RowTable(Any[r for r in rows], colindex)
"""
    RowTable(rows::Vector, colindex::CIndex)

Construct a `RowTable` from `rows` and a list of column names `colindex`.
"""
RowTable(rows::Vector, colindex::CIndex) = RowTable(Any[rows...], colindex)

### Access

@inline cindex(rt::RowTable) = rt.colindex
@inline colindex(rt::RowTable, ci) = cindex(rt)[ci]  # use this, or not

"""
    rows(rt::RowTable)

Return the rows of `rt` as a Vector. `rt` itself is iteratble and
iterates over rows. But, iterating over `rows(rt)` is faster.
"""
@inline rows(rt::RowTable) = rt.rows

# Annotating the method above with return type Vector makes it 1000x slower in v0.6,
# and only a bit less slower in v0.7.
# @inline rows(rt::RowTable)::Vector = rt.rows

"""
    names(rt::RowTable)

Return a list of the column names of `rt`.
"""
Base.names(rt::RowTable) = _names(cindex(rt))
_names(rt::RowTable) = _names(cindex(rt))

"""
    size(rt::RowTable)

Return the number of rows and the number columns in `rt`.
"""
Base.size(rt::RowTable) = _numberofrows(rt), _numberofcols(rt)
_numberofcols(rt::RowTable) = length(_names(rt))
_numberofrows(rt::RowTable) = length(rows(rt))

## TODO: make sure this is optimized if n is known at compile time
## (And even if not known)
Base.size(rt::RowTable, n::Integer) = n == 1 ? _numberofrows(rt) : n == 2 ? _numberofcols(rt) :
    error(ArgumentError, ": RowTables have only two dimensions")

### Equality

Base.:(==)(rt1::RowTable, rt2::RowTable) = (cindex(rt1) == cindex(rt2) && rows(rt1) == rows(rt2))

### Constructors

newrows(n::Integer=0) = Vector{Any}(undef, n)

## We did not use the type information afterall.
const _NameTypes = Union{AbstractVector{S} where S<:Union{Symbol, AbstractString}, Tuple}

function _RowTable(a::AbstractVector{T}, keynames; tuples::Bool=false) where {T <: AbstractDict}
    isempty(a) && return RowTable(newrows(), CIndex(map(Symbol, keynames))) # JSON keys are strings
    l = length(first(a))
    all(x -> length(x) == l, a) || throw(DimensionMismatch("All dictionaries must be of the same length"))
    if tuples
        RowTable([dict_to_named_tuple(a[i], keynames) for i in LinearIndices(a)], CIndex(map(Symbol, keynames)))
    else
        RowTable([map(x -> a[i][x], keynames) for i in LinearIndices(a)], CIndex(map(Symbol, keynames)))
    end
end

# FIXME: ensue symbols as keys
function dict_to_named_tuple(d::AbstractDict, keynames)
    datatypes = Tuple{[typeof(d[k]) for k in keynames]...,}
    return NamedTuple{(keynames...,), datatypes}(([d[k] for k in keynames]...,))
end

# Construct from Array of Dicts
RowTable(a::AbstractVector{T}, keynames::_NameTypes; tuples=false) where {T<:AbstractDict} = _RowTable(a, keynames, tuples=tuples)

function _RowTable(::Type{T}, a::AbstractVector, keynames) where T <: AbstractArray
    all(x -> isa(x, AbstractArray), a) || error("Not all elements are arrays")  # They don't have to be. Just not dicts
    RowTable(a, CIndex(keynames))
end

function _RowTable(::Type{T}, a::AbstractVector, keynames) where T <: Tuple
    RowTable(a, CIndex(keynames))
end

# v0.7 requires collect (or something else) here to avoid constructing a Set, which prevents indexing
_RowTable(::Type{T}, a::AbstractVector) where T <: AbstractDict  = _RowTable(T, a, collect(keys(first(a))))

function _RowTable(::Type{T}, a::AbstractVector, keynames; tuples=false) where T <: AbstractDict
    all(x -> isa(x, AbstractDict), a) || error("Not all elements are dictionaries")
    _RowTable(a, keynames, tuples=tuples)
end

function RowTable(a::AbstractVector)
    isempty(a) && return RowTable()
    _RowTable(typeof(first(a)), a)
end

"""
    RowTable(a::AbstractVector, keynames)

Construct a `RowTable` from a vector of rows `a`, and `keynames`.
If `keynames` is a `Vector`, then the rows are `Vectors`.
If `keynames` is a `Tuple`, they are `Tuple`s.

`RowTable([], keynames)` constructs a `RowTable` with no rows.
"""
function RowTable(a::AbstractVector, keynames::_NameTypes)
    isempty(a) && return RowTable([], CIndex(keynames))
    _RowTable(typeof(first(a)), a, keynames)
end

"""
    RowTable(cols=[], names=[], tuples=false)

Return a `RowTable` constructed from the collection of columns `cols`,
with `names`. If `tuples=true`, the rows are named tuples. Otherwise,
they are `Array`s.
"""
function RowTable(; cols=[], names=[], tuples=false)
    if isempty(cols)
        isempty(names) && return RowTable(newrows(), CIndex())
        return RowTable(newrows(), CIndex(names))
    end
    length(cols) == length(names) || throw(DimensionMismatch("Number of columns and number of names must be equal."))
    isempty(cols[1]) && throw(DimensionMismatch("Support for empty columns not implmented."))
    columnlength = length(cols[1])
    for c in cols
        (length(c) == columnlength) || throw(DimensionMismatch("Columns must have the same length."))
    end
    rows = tuples ? namedtuple_rows(cols, names, columnlength) : array_rows(cols, columnlength)
    return RowTable(rows, CIndex(names))
end

function namedtuple_rows(cols, names, columnlength)
    datatypetup = (collect(c[1] for c in cols)...,) |> typeof
    nametup = (names...,)
    namedtupletype = NamedTuple{nametup}{datatypetup}
    rows = newrows(columnlength)
    for i in 1:columnlength
        rows[i] = namedtupletype(([c[i] for c in cols]...,))
    end
    return rows
end

function array_rows(cols, columnlength)
    rows = newrows(columnlength)
    for i in 1:columnlength
        rows[i] = [c[i] for c in cols]
    end
    return rows
end

"""
    RowTable(df::DataFrames.DataFrame; tuples=false)

Construct a `RowTable` from a `DataFrame`. If `tuples=true`
then the constructed rows are `Tuples` rather than `Vectors`.
"""
function RowTable(df::DataFrames.DataFrame; tuples=false)
    (nr, nc) = size(df)
    arr = newrows()
    if tuples
        @inbounds for ri in 1:nr
            push!(arr, ([df[ri, ci] for ci in 1:nc]...,))
        end
    else
        @inbounds for ri in 1:nr
            push!(arr, [df[ri, ci] for ci in 1:nc])
        end
    end
    RowTable(arr, copy(names(df)))
end

##############################################################################
##
## getindex()
##
##############################################################################

## This allows single methods to handle both Int and Symbol indices.
## But, it allows repeated mapping of the same Symbol, which is inefficient.
## So, this might not be used much
const ColInd = Union{Integer, Symbol}

## A single index is a interpreted as a column index, consistent with DataFrames
Base.getindex(rt::RowTable, cinds) = rt[:,cinds]

## Return element in a single cell
Base.getindex(rt::RowTable, ri::Integer, ci::Symbol) = rows(rt)[ri][cindex(rt)[ci]]

# If above is called in a loop with symbol arg, using below is faster

Base.getindex(rt::RowTable, ri::Integer, ci::Integer) = rows(rt)[ri][ci]

## Return a slice of a column as a Vector
function Base.getindex(rt::RowTable, ri::AbstractVector, ci::ColInd)
    ind = colindex(rt, ci) # do this so symbol mapping is only done once
    [rt[i, ind] for i in ri]
end

## Return a slice of a row as a Vector
Base.getindex(rt::RowTable, ri::Integer, cis::AbstractVector{T}) where {T <: Symbol} =
    rows(rt)[ri][[cindex(rt)[ci] for ci in cis]]

Base.getindex(rt::RowTable, ri::Integer, cis::AbstractVector{T}) where {T <: Integer} =
    rows(rt)[ri][cis]

### Return a row, vector or Tuple
Base.getindex(rt::RowTable, ri::Integer, ::Colon) = rt.rows[ri]

## Return slice as RowTable
## Following method calls the next method with integer arguments
Base.getindex(rt::RowTable, ri::AbstractVector{T}, ci::AbstractVector{V}) where {T<:Integer, V<:Symbol} =
    Base.getindex(rt, ri, [cindex(rt)[s] for s in ci])

## Return rectangular slice in both dimensions as RowTable
function Base.getindex(rt::RowTable, ri::AbstractVector, ci::AbstractVector{T}) where T<:Integer
    ar = newrows(length(ri))
    for (i, ind) in enumerate(ri)
        ar[i] = rows(rt)[ind][ci]
    end
    RowTable(ar, CIndex(cindex(rt).names[ci]))
end

Base.getindex(rt::RowTable, ::Colon, ci) = rt[1:length(rows(rt)), ci]

Base.getindex(rt::RowTable, ri::AbstractVector, ::Colon) = RowTable(rows(rt)[ri], rt.colindex)

##############################################################################
##
## setindex!()
##
##############################################################################

## Set a single element
Base.setindex!(rt::RowTable, val, ri::Integer, ci::Integer) = (rows(rt)[ri][ci] = val)
Base.setindex!(rt::RowTable, val, ri::Integer, ci::Symbol) = (rows(rt)[ri][cindex(rt)[ci]] = val)

### Convert

"""
    columns(rt::RowTable)::Vector

Return the columns of `rt`.
"""
function columns(rt::RowTable)
    (nr, nc) = size(rt)
    @inbounds colarr =  [newrows(nr) for i in 1:nc] # misusing newrows for columns here
    return _columns!(rt, colarr)
end

function columnstyped(rt::RowTable)
    (nr, nc) = size(rt)
    @inbounds colarr =  [Vector{typeof(rt[1, i])}(undef, nr) for i in 1:nc]
    return _columns!(rt, colarr)
end

@inline function _columns!(rt, colarr)
    (nr, nc) = size(rt)
    for rowind in 1:nr
      @inbounds row = rows(rt)[rowind]
        for colind in 1:nc
        @inbounds colarr[colind][rowind] = row[colind]
        end
    end
    return colarr
end

"""
    DataFrame(rt::RowTable; typed=false)

Convert `rt` to a `DataFrame`. If `typed` is `true`,
then the eltype of each column is the type of the first
element in the column. This will raise an error if the
elements are in fact not of the same type.
"""
function DataFrames.DataFrame(rt::RowTable; typed=false)
    cols = (typed ? columnstyped(rt) : columns(rt))
    DataFrames.DataFrame(cols,_names(rt))
end

### Copy

Base.copy(rt::RowTable) = RowTable(copy(rows(rt)), copy(cindex(rt)))
Base.deepcopy(rt::RowTable) = RowTable(deepcopy(rows(rt)), deepcopy(cindex(rt)))

### Iterate over rows
## DataFrames does not define iterating over a DataFrame, rather requires
## specifying rows or columns

for f in (:length, :iterate)
    @eval begin
        (Base.$f)(rt::RowTable) = (Base.$f)(rows(rt))
    end
end

for f in (:iterate,)
    @eval begin
        (Base.$f)(rt::RowTable, args...) = (Base.$f)(rows(rt), args...)
    end
end

### Transform

for f in (:deleteat!, :push!, :pushfirst!, :insert!, :append!, :prepend!, :splice!, :permute!)
    @eval begin
        (Base.$f)(rt::RowTable, args...) = (($f)(rows(rt), args...); rt)
    end
end

for f in (:pop!, :popfirst!)
    @eval begin
        (Base.$f)(rt::RowTable, args...) = ($f)(rows(rt), args...)
    end
end

permute!(rt::RowTable, p::AbstractVector) = (permute!(rows(rt), p); rt)
permute(rt::RowTable, p::AbstractVector) = permute!(copy(rt), p)

Random.shuffle!(rng::Random.AbstractRNG, rt::RowTable) = (shuffle!(rng, rows(rt)); rt)
Random.shuffle!(rt::RowTable) = (shuffle!(rows(rt)); rt)
Random.shuffle(rt::RowTable) = shuffle!(copy(rt))
Random.shuffle(rng::Random.AbstractRNG, rt::RowTable) = shuffle!(rng, copy(rt))

DataFrames.rename!(rt::RowTable, d) = (rename!(cindex(rt), d); rt)
DataFrames.rename(rt::RowTable, d) = rename!(copy(rt), d)


text:
@testset "Locally free class group of group algebras" begin
  G = small_group(8, 4)
  A = AlgAss(AlgGrp(FlintQQ, G))[1]
  O = Order(A, basis(A))
  C = Hecke.locally_free_class_group(O)
  @test C.snf == fmpz[ 2 ]

  G = small_group(10, 2)
  A = AlgAss(AlgGrp(FlintQQ, G))[1]
  O = Order(A, basis(A))
  C = Hecke.locally_free_class_group(O)
  @test C.snf == fmpz[]

  G = small_group(12, 3)
  A = AlgAss(AlgGrp(FlintQQ, G))[1]
  O = Order(A, basis(A))
  C = Hecke.locally_free_class_group(O)
  @test C.snf == fmpz[]
end

@testset "Locally free class group of matrix algebras" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 47
  K, a = number_field(f, "a")
  A = AlgAss(MatrixAlgebra(K, 2))
  A, _, _ = Hecke.restrict_scalars(A, FlintQQ)
  O = MaximalOrder(A)
  C = Hecke.locally_free_class_group(O) # == class_group(K)
  @test C.snf == fmpz[ 5 ]

  f = x^2 + 26
  K, a = number_field(f, "a")
  A = AlgAss(MatrixAlgebra(K, 2))
  A, _, _ = Hecke.restrict_scalars(A, FlintQQ)
  O = MaximalOrder(A)
  C = Hecke.locally_free_class_group(O) # == class_group(K)
  @test C.snf == fmpz[ 6 ]
end


text:
immutable Downhill{G<:Function} <: Policy
    gradient::G
end
Downhill(p::PowseekerProblem) = Downhill(mdp(p).gradient)

function action(p::Downhill, s::SkierState)
    grad = p.gradient(s.pos)
    steepest = atan2(-grad[2], -grad[1])
    return steepest - s.psi
end

immutable Topout{G<:Function} <: Policy
    mdp::PowseekerMDP{G}
    tol::Float64
end
Topout(p::PowseekerProblem, tol::Float64) = Topout(mdp(p), tol)

function action(p::Topout, s::SkierState)
    grad = p.mdp.gradient(s.pos)
    if s.time >= p.mdp.duration || norm(grad) < p.tol # topped out
        return action(Downhill(p.mdp.gradient), s)
    else
        up = atan2(grad[2], grad[1])
        return up-s.psi
    end
end

immutable SkiOver{G<:Function} <: Policy
    mdp::PowseekerMDP{G}
end
SkiOver(p::PowseekerProblem) = SkiOver(mdp(p))

function action(p::SkiOver, s::SkierState)
    grad = p.mdp.gradient(s.pos)
    if s.time >= p.mdp.duration
        return action(Downhill(p.mdp.gradient), s)
    elseif dot(grad, [cos(s.psi), sin(s.psi)]) > 0 # facing uphill
        up = atan2(grad[2], grad[1])
        return up-s.psi
    else
        down = atan2(-grad[2], -grad[1])
        return down-s.psi
    end
end

immutable RandomlyCheckGPS{P<:Policy, RNG<:AbstractRNG} <: Policy
    p::P
    p_check::Float64
    rng::RNG
end

function action(p::RandomlyCheckGPS, b)
    check = rand(p.rng) < p.p_check
    if check
        return GPSOrAngle(true, 0.0)
    else
        s = rand(p.rng, b)
        return GPSOrAngle(false, action(p.p, s))
    end
end

#=
@with_kw immutable Center <: Policy
    rng::MersenneTwister = MersenneTwister(14)
end

function action(::Center, s::SkierState)
    center = atan2
end

function action(p::Center, c::ParticleCollection)
    s = rand(p.rng, c)
    return GPSOrAngle(false, action(p, s))
end
=#


text:
export ConstantOperator, IdentityOperator, BasisFunctional

struct ConstantOperator{T,DS} <: Operator{T}
    λ::T
    space::DS
    ConstantOperator{T,DS}(c::Number,sp::DS) where {T,DS} = new{T,DS}(convert(T,c),sp)
    ConstantOperator{T,DS}(L::UniformScaling,sp::DS) where {T,DS} = new{T,DS}(convert(T,L.λ),sp)
end

ConstantOperator(::Type{T},c,sp::Space) where {T} = ConstantOperator{T,typeof(sp)}(c,sp)
ConstantOperator(::Type{T},c) where {T} = ConstantOperator(T,c,UnsetSpace())
ConstantOperator(c::Number,sp::Space) = ConstantOperator(typeof(c),c,sp)
ConstantOperator(c::Number) = ConstantOperator(typeof(c),c)
ConstantOperator(L::UniformScaling) = ConstantOperator(L.λ)
ConstantOperator(L::UniformScaling,sp::Space) = ConstantOperator(L.λ,sp)
IdentityOperator() = ConstantOperator(1.0)
IdentityOperator(S::Space) = ConstantOperator(1.0,S)

for OP in (:domainspace,:rangespace)
    @eval $OP(C::ConstantOperator) = C.space
end

promotedomainspace(C::ConstantOperator,sp::Space) = ConstantOperator(C.λ,sp)
promoterangespace(C::ConstantOperator,sp::Space,cursp::UnsetSpace) = ConstantOperator(C.λ,sp)

bandwidths(T::ConstantOperator) = 0,0

isbandedblockbanded(::ConstantOperator) = true
blockbandwidths(::ConstantOperator) = 0,0
subblockbandwidths(::ConstantOperator) = 0,0

getindex(C::ConstantOperator,k::Integer,j::Integer) =
    k==j ? eltype(C)(C.λ) : zero(eltype(C))

==(C1::ConstantOperator, C2::ConstantOperator) = C1.λ==C2.λ

function convert(::Type{Operator{T}}, C::ConstantOperator) where T
    if T == eltype(C)
        C
    else
        ConstantOperator{T,typeof(C.space)}(C.λ,C.space)
    end
end

# zero needs to be different since it can take a space to
# a ConstantSpace, in creating functionals
convert(::Type{Operator{T}}, x::Number) where {T} =
    x==0 ? ZeroOperator(T) : Multiplication(convert(T,x))
convert(::Type{Operator{T}}, L::UniformScaling) where {T} =
    ConstantOperator(T,L.λ)

convert(::Type{Operator},n::Number) = convert(Operator{typeof(n)}, n)
convert(::Type{Operator},L::UniformScaling) = ConstantOperator(L.λ)

## Algebra

for op in (:+,:-,:*)
    @eval ($op)(A::ConstantOperator,B::ConstantOperator) = ConstantOperator($op(A.λ,B.λ))
end

## Basis Functional

struct BasisFunctional{T,DS} <: Operator{T}
    k::Int
    space::DS
end

@functional BasisFunctional

BasisFunctional(k,sp) = BasisFunctional{Float64,typeof(sp)}(k, sp)
BasisFunctional(k) = BasisFunctional(k, ℓ⁰)

bandwidths(B::BasisFunctional) = 0,B.k-1
domainspace(B::BasisFunctional) = B.space

convert(::Type{Operator{T}},B::BasisFunctional) where {T} = BasisFunctional{T,typeof(B.space)}(B.k,B.space)

Base.getindex(op::BasisFunctional{T},k::Integer) where {T} = (k==op.k) ? one(T) : zero(T)
Base.getindex(op::BasisFunctional{T},k::AbstractRange) where {T} = convert(Vector{T},k.==op.k)

struct FillFunctional{T} <: Operator{T}
    λ::T
end

@functional FillFunctional

domainspace(B::FillFunctional) = ℓ⁰

Base.getindex(op::FillFunctional,k::Integer)=op.λ
Base.getindex(op::FillFunctional,k::AbstractRange)=fill(op.λ,length(k))

## Zero is a special operator: it makes sense on all spaces, and between all spaces

struct ZeroOperator{T,S,V} <: Operator{T}
    domainspace::S
    rangespace::V
end

ZeroOperator(::Type{T},d::Space,v::Space) where {T} = ZeroOperator{T,typeof(d),typeof(v)}(d,v)
ZeroOperator(::Type{T},S::Space) where {T} = ZeroOperator(T,S,ZeroSpace(S))
ZeroOperator(d::Space,v::Space) = ZeroOperator(Float64,d,v)
ZeroOperator(S::Space) = ZeroOperator(S,ZeroSpace(S))
ZeroOperator() = ZeroOperator(UnsetSpace(),UnsetSpace())
ZeroOperator(::Type{T}) where {T} = ZeroOperator(T,UnsetSpace(),UnsetSpace())

convert(::Type{Operator{T}},Z::ZeroOperator) where {T} =
    ZeroOperator(T,Z.domainspace,Z.rangespace)

Base.zero(A::Operator) = ZeroOperator(domainspace(A),rangespace(A))

domainspace(Z::ZeroOperator)=Z.domainspace
rangespace(Z::ZeroOperator)=Z.rangespace

bandwidths(T::ZeroOperator) = -720,-720   # 6!
blockbandwidths(T::ZeroOperator) = -720,-720   # 6!
subblockbandwidths(T::ZeroOperator) =  -720, -720

isbandedblockbandedabove(::ZeroOperator) = true
isbandedblockbandedbelow(::ZeroOperator) = true

getindex(C::ZeroOperator,k::Integer) = zero(eltype(C))
getindex(C::ZeroOperator,k::Integer,j::Integer) = zero(eltype(C))

promotedomainspace(Z::ZeroOperator,sp::UnsetSpace) = Z
promoterangespace(Z::ZeroOperator,sp::UnsetSpace) = Z
promotedomainspace(Z::ZeroOperator,sp::Space) = ZeroOperator(sp,rangespace(Z))
promoterangespace(Z::ZeroOperator,sp::Space) = ZeroOperator(domainspace(Z),sp)

isconstop(::Union{ZeroOperator,ConstantOperator}) = true
isconstop(S::SpaceOperator) = isconstop(S.op)
isconstop(_) = false

iszeroop(::ZeroOperator) = true
iszeroop(A::ConstantOperator) = A.λ==0.0
iszeroop(A) = false

convert(::Type{T},::ZeroOperator) where {T<:Number} = zero(T)
convert(::Type{T},C::ConstantOperator) where {T<:Number} = convert(T,C.λ)
convert(::Type{T},S::SpaceOperator) where {T<:Number} = convert(T,S.op)

## SubOperator convert
## Special case for ZeroOperator
for TYP in (:RaggedMatrix,:Matrix,:BandedMatrix,
            :BlockBandedMatrix,:BandedBlockBandedMatrix)
    @eval begin
        $TYP(S::SubOperator{T,ZO}) where {T,ZO<:ZeroOperator} = $TYP(Zeros, S)
        $TYP(S::SubOperator{T,ZO,NTuple{2,UnitRange{Int}}}) where {T,ZO<:ZeroOperator} = $TYP(Zeros, S)
    end
end


text:
using ExecutableSpecifications:
    FromMacroStepDefinitionMatcher, executefeature, Executor, CompositeStepDefinitionMatcher,
    ColorConsolePresenter, present, ResultAccumulator, accumulateresult, issuccess, featureresults
using ExecutableSpecifications.Gherkin: parsefeature, Given, When, Then, Feature

function allfileswithext(path::String, extension::String)
    [filename for filename in readdir(path)
              if isfile(joinpath(path, filename)) &&
                 splitext(joinpath(path, filename))[2] == extension]
end

function runspec()
    # Find all step definition files and all feature files.
    stepfiles = allfileswithext("features/steps", ".jl")
    featurefiles = allfileswithext("features", ".feature")

    # Read all step definition files.
    matchers = FromMacroStepDefinitionMatcher[]
    for filename in stepfiles
        fullpath = joinpath("features/steps", filename)
        push!(matchers, FromMacroStepDefinitionMatcher(readstring(fullpath); filename=fullpath))
    end
    matcher = CompositeStepDefinitionMatcher(matchers...)

    # Read all feature files.
    features = Feature[]
    for filename in featurefiles
        featureresult = parsefeature(readstring(joinpath("features", filename)))
        push!(features, featureresult.value)
    end

    # The presenter will print all scenario steps as they are executed.
    presenter = ColorConsolePresenter()
    executor = Executor(matcher, presenter)
    accumulator = ResultAccumulator()

    for feature in features
        result = executefeature(executor, feature)
        accumulateresult(accumulator, result)
    end

    println()

    #
    # Present number of scenarios that succeeded and failed for each feature
    #
    results = featureresults(accumulator)

    # Find the longest feature name, so we can align the result table.
    maxfeature = maximum(length(r.feature.header.description) for r in results)

    featureprefix = "  Feature: "
    print_with_color(:white, " " ^ (length(featureprefix) + maxfeature + 1), "| Success | Failure\n")
    for r in results
        linecolor = r.n_failure == 0 ? :green : :red
        print_with_color(linecolor, featureprefix, rpad(r.feature.header.description, maxfeature))
        print_with_color(:white, " | ")
        print_with_color(:green, rpad("$(r.n_success)", 7))
        print_with_color(:white, " | ")
        print_with_color(linecolor, rpad("$(r.n_failure)", 7), "\n")
    end

    println()

    istotalsuccess = issuccess(accumulator)
    if istotalsuccess
        println("SUCCESS")
    else
        println("FAILURE")
    end

    istotalsuccess
end


text:
function jlmie_Enf(nmat,radius,lbd0,nenv,x,y,z,nmax::Int=-1)
    if nmax == -1
        nmax = largenmax(jlmie_mx(nmat,radius,lbd0,nenv)[2])
    end

end

function jlmie_Enf1_n(nmat,radius,lbd0,nenv,r,theta,phi,n::Int=1)
# implementation of eq.(4.40) of B&H
    m = nmat./nenv
    lbd = lbd0./nenv  # effective wavelength in the environment
    lbd1 = lbd0./nmat  # effective wavelength in the sphere
    k = 2*pi./lbd  # free space wavenumber in the environment
    k1 = 2*pi./lbd1  # free space wavenumber in the sphere
    x = k.*radius
    x1 = k1.*radius
    mu = cos.(theta)
    _,_,cn,dn = jlmie_abcd(m,x,n)
    pin,taun = jlmie_pt(cos.(theta),n)

    jx1 = sbesselj.(n,x1)
    d1x1jx1 = x1.*sbesselj.(n-1,x1) .- n.*sbesselj.(n,x1)
    # d1mxjmx = m.*x.*sbesselj.(n-1,m.*x) .- n.*sbesselj.(n,m.*x)

    E0 = 1.0;
    En = (1im)^n *E0*(2*n+1)/(n*(n+1))

    M1o1n_r = 0
    M1o1n_t = cos.(phi).*pin.*jx1
    M1o1n_p = -1*sin.(phi).*taun.*jx1
    # M1e1n_r = 0
    # M1e1n_t = -1*sin.(phi).*pin.*jx1
    # M1e1n_p = -1*cos.(phi).*taun.*jx1
    # N1o1n_r = sin.(phi).*n.*(n+1).*sin.(theta).*pin.*jx1./x1
    # N1o1n_t = sin.(phi).*tau.*d1x1jx1./x1
    # N1o1n_p = cos.(phi).*pin.*d1x1jx1./x1
    N1e1n_r = cos.(phi).*n.*(n+1).*sin.(theta).*pin.*jx1./x1
    N1e1n_t = cos.(phi).*tau.*d1x1jx1./x1
    N1e1n_p = -1*sin.(phi).*pin.*d1x1jx1./x1

    E1r_n = En.*(cn.*M1o1n_r .- 1im.*dn.*N1e1n_r)
    E1t_n = En.*(cn.*M1o1n_t .- 1im.*dn.*N1e1n_t)
    E1p_n = En.*(cn.*M1o1n_p .- 1im.*dn.*N1e1n_p)
    return E1r_n,E1t_n,E1p_n
end

function jlmie_Enfs_n(nmat,radius,lbd0,nenv,r,theta,phi,n::Int=1)
    # implementation of eq.(4.45) of B&H
    m = nmat./nenv
    lbd = lbd0./nenv  # effective wavelength in the environment
    k = 2*pi./lbd  # free space wavenumber in the environment
    x = k.*radius
    mu = cos.(theta)
    an,bn,_,_ = jlmie_abcd(m,x,n)
    pin,taun = jlmie_pt(cos.(theta),n)

    hx = shankelh1.(n,x)
    d1xhx = x.*shankelh1.(n-1,x) .- n.*shankelh1.(n,x)

    E0 = 1.0;
    En = (1im)^n *E0*(2*n+1)/(n*(n+1))

    M3o1n_r = 0
    M3o1n_t = cos.(phi).*pin.*hx
    M3o1n_p = -1*sin.(phi).*taun.*hx
    # M3e1n_r = 0
    # M3e1n_t = -1*sin.(phi).*pin.*hx
    # M3e1n_p = -1*cos.(phi).*taun.*hx
    # N3o1n_r = sin.(phi).*n.*(n+1).*sin.(theta).*pin.*hx./x
    # N3o1n_t = sin.(phi).*tau.*d1xhx./x
    # N3o1n_p = cos.(phi).*pin.*d1xhx./x
    N3e1n_r = cos.(phi).*n.*(n+1).*sin.(theta).*pin.*hx./x
    N3e1n_t = cos.(phi).*tau.*d1xhx./x
    N3e1n_p = -1*sin.(phi).*pin.*d1xhx./x

    Esr_n = En.*(1im.*an.*N3e1n_r .- bn.*M3o1n_r)
    Est_n = En.*(1im.*an.*N3e1n_t .- bn.*M3o1n_t)
    Esp_n = En.*(1im.*an.*N3e1n_p .- bn.*M3o1n_p)
    return Esr_n,Est_n,Esp_n
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
const stdbravais = Dict{String, LatticeParameter}()
const stdunitcells = Dict{String, LatticeParameter}()
const stdlattices = Dict{String, LatticeParameter}()

## 1D bravais
stdbravais["chain"] = P(:name =>"chain",
                        :dimension => 1,
                        :parameters => [(:a,1.0)],
                        :basis => :(reshape([a],1,1)))

## 2D bravais
stdbravais["orthorhombic2d"] = P(:name => "orthorhombic2d",
                                 :dimension => 2,
                                 :parameters => [(:a,1.0), (:b,1.0)],
                                 :basis => :([a 0
                                              0 b]))
stdbravais["hexagonal2d"] = P(:name => "hexagonal2d",
                              :dimension => 2,
                              :parameters => [(:a,1.0)],
                              :basis => :([a  -a/2
                                           0  a*sqrt(3.0)/2]))
stdbravais["monoclinic2d"] = P(:name => "monoclinic2d",
                               :dimension => 2,
                               :parameters => [(:a,1.0), (:b,1.0), (:theta, 90.0)],
                               :basis => :([a b*cosd(theta)
                                            0 b*sind(theta)]))

## 3D bravais
stdbravais["orthorhombic3d"] = P(:name => "orthorhombic3d",
                                 :dimension => 3,
                                 :parameters => [(:a,1.0),(:b,1.0),(:c,1.0)],
                                 :basis => :([a 0 0
                                              0 b 0
                                              0 0 c]))

## 1D unitcells
stdunitcells["simple1d"] = P(:name => "simple1d",
                             :dimension => 1,
                             :sites => [P(:id=>1,
                                          :sitetype=>1,
                                          :coord=>[0.0]
                                         )],
                             :bonds => [P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0]),
                                          :target=>P(:id=>1, :offset=>[1]),
                                         )])

stdunitcells["bond-alternating simple1d"] = P(:name => "bond-alternating simple1d",
                                              :dimension => 1,
                                              :sites => [P(:id=>1,
                                                           :sitetype=>1,
                                                           :coord=>[0.0],
                                                          ),
                                                         P(:id=>2,
                                                           :sitetype=>2,
                                                           :coord=>[0.5],
                                                          ),
                                                        ],
                                              :bonds => [P(:bondtype=>1,
                                                           :source=>P(:id=>1, :offset=>[0]),
                                                           :target=>P(:id=>2, :offset=>[0]),
                                                          ),
                                                         P(:bondtype=>2,
                                                           :source=>P(:id=>2, :offset=>[0]),
                                                           :target=>P(:id=>1, :offset=>[1]),
                                                           )])

## 2D unitcells
stdunitcells["simple2d"] = P(:name => "simple2d",
                             :dimension => 2,
                             :sites => [P(:id=>1,
                                          :sitetype=>1,
                                          :coord=>[0.0, 0.0],
                                         )],
                             :bonds => [P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0,0]),
                                          :target=>P(:id=>1, :offset=>[1,0]),),
                                        P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0,0]),
                                          :target=>P(:id=>1, :offset=>[0,1]),),
                                       ])
stdunitcells["triangular cell"] = P(:name => "triangular cell",
                                    :dimension => 2,
                                    :sites => [P(:id=>1,
                                                 :sitetype=>1,
                                                 :coord=>[0.0,0.0],
                                                )],
                                    :bonds => [P(:bondtype=>1,
                                                 :source=>P(:id=>1, :offset=>[0,0]),
                                                 :target=>P(:id=>1, :offset=>[1,0]),),
                                               P(:bondtype=>1,
                                                 :source=>P(:id=>1, :offset=>[0,0]),
                                                 :target=>P(:id=>1, :offset=>[0,1]),),
                                               P(:bondtype=>1,
                                                 :source=>P(:id=>1, :offset=>[0,0]),
                                                 :target=>P(:id=>1, :offset=>[1,1]),),
                                              ])
stdunitcells["honeycomb cell"] = P(:name => "honeycomb cell",
                                   :dimension => 2,
                                   :sites => [P(:id=>1,
                                                :sitetype=>1,
                                                :coord=>[0.0, 0.0],
                                               ),
                                              P(:id=>2,
                                                :sitetype=>2,
                                                :coord=> (2//3).*[0.5, 1.0],
                                               )],
                                   :bonds => [P(:bondtype=>1,
                                                :source=>P(:id=>1, :offset=>[0,0]),
                                                :target=>P(:id=>2, :offset=>[0,0]),),
                                              P(:bondtype=>2,
                                                :source=>P(:id=>2, :offset=>[0,0]),
                                                :target=>P(:id=>1, :offset=>[1,1]),),
                                              P(:bondtype=>3,
                                                :source=>P(:id=>2, :offset=>[0,0]),
                                                :target=>P(:id=>1, :offset=>[0,1]),),
                                             ])

## 3D unitcells
stdunitcells["simple3d"] = P(:name => "simple3d",
                             :dimension => 3,
                             :sites => [P(:id=>1,
                                          :sitetype=>1,
                                          :coord=>[0.0, 0.0, 0.0],
                                         )],
                             :bonds => [P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0,0,0]),
                                          :target=>P(:id=>1, :offset=>[1,0,0]),),
                                        P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0,0,0]),
                                          :target=>P(:id=>1, :offset=>[0,1,0]),),
                                        P(:bondtype=>1,
                                          :source=>P(:id=>1, :offset=>[0,0,0]),
                                          :target=>P(:id=>1, :offset=>[0,0,1]),),
                                       ])

## 1D lattices
stdlattices["chain lattice"] = P(:name => "chain lattice",
                                 :dimension => 1,
                                 :bravais => "chain",
                                 :unitcell => "simple1d",
                                 :parameters => [],
                                 :periodic => [true],
                                )

stdlattices["bond-alternating chain lattice"] = P(:name => "bond-alternating chain lattice",
                                 :dimension => 1,
                                 :bravais => "chain",
                                 :unitcell => "bond-alternating simple1d",
                                 :parameters => [(:a, 2.0)],
                                 :periodic => [true],
                                )

## 2D lattices
stdlattices["square lattice"] = P(:name => "square lattice",
                                  :dimension => 2,
                                  :bravais => "orthorhombic2d",
                                  :unitcell => "simple2d",
                                  :parameters => [],
                                  :periodic => [true, true],
                                 )
stdlattices["triangular lattice"] = P(:name => "triangular lattice",
                                      :dimension => 2,
                                      :bravais => "hexagonal2d",
                                      :unitcell => "triangular cell",
                                      :parameters => [],
                                      :periodic => [true, true],
                                     )
stdlattices["honeycomb lattice"] = P(:name => "honeycomb lattice",
                                     :dimension => 2,
                                     :bravais => "hexagonal2d",
                                     :unitcell => "honeycomb cell",
                                     :parameters => [(:a, sqrt(3.0))],
                                     :periodic => [true, true],
                                    )
stdlattices["ladder"] = P(:name => "ladder",
                          :dimension => 2,
                          :bravais => "orthorhombic2d",
                          :unitcell => "simple2d",
                          :parameters => [],
                          :periodic => [true, false],
                         )

## 3D lattices
stdlattices["cubic lattice"] = P(:name => "cubic lattice",
                                 :dimension => 3,
                                 :bravais => "orthorhombic3d",
                                 :unitcell => "simple3d",
                                 :periodic => [true, true, true],
                                )

## special lattices
stdlattices["fully connected graph"] = P(:name => "fully connected graph",
                                           :dim => 1,
                                           :generator => generate_fully_connected_graph,
                                          )


text:
### Base HDF5 structs ###

# High-level reference handler
struct Reference
  r::API.hobj_ref_t
end
Reference() = Reference(API.HOBJ_REF_T_NULL) # NULL reference to compare to
Base.cconvert(::Type{Ptr{T}}, ref::Reference) where {T<:Union{Reference,API.hobj_ref_t,Cvoid}} = Ref(ref)
Base.:(==)(a::Reference, b::Reference) = a.r == b.r
Base.hash(x::Reference, h::UInt) = hash(x.r, h)

function Reference(parent::Union{File,Group,Dataset}, name::AbstractString)
  ref = Ref{API.hobj_ref_t}()
  API.h5r_create(ref, checkvalid(parent), name, API.H5R_OBJECT, -1)
  return Reference(ref[])
end

# Single character types
# These are needed to safely handle VLEN objects
abstract type CharType <: AbstractString end

struct ASCIIChar <: CharType
  c::UInt8
end
Base.length(c::ASCIIChar) = 1

struct UTF8Char <: CharType
  c::UInt8
end
Base.length(c::UTF8Char) = 1

chartype(::Type{String}) = ASCIIChar
stringtype(::Type{ASCIIChar}) = String
stringtype(::Type{UTF8Char}) = String

cset(::Type{<:AbstractString}) = API.H5T_CSET_UTF8
cset(::Type{UTF8Char}) = API.H5T_CSET_UTF8
cset(::Type{ASCIIChar}) = API.H5T_CSET_ASCII

const BitsType = Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float32,Float64}
const ScalarType = Union{BitsType,Reference}

# VLEN objects
struct VLen{T}
  data::Array
end
VLen(strs::Array{S}) where {S<:String} = VLen{chartype(S)}(strs)
VLen(A::Array{Array{T}}) where {T<:ScalarType} = VLen{T}(A)
VLen(A::Array{Array{T,N}}) where {T<:ScalarType,N} = VLen{T}(A)
function Base.cconvert(::Type{Ptr{Cvoid}}, v::VLen)
  len = length(v.data)
  h = Vector{API.hvl_t}(undef, len)
  for ii in 1:len
      d = v.data[ii]
      p = unsafe_convert(Ptr{UInt8}, d)
      h[ii] = API.hvl_t(length(d), p)
  end
  return h
end

datatype(A::VLen{T}) where {T<:ScalarType} = Datatype(API.h5t_vlen_create(hdf5_type_id(T)))
function datatype(str::VLen{C}) where {C<:CharType}
    type_id = API.h5t_copy(hdf5_type_id(C))
    API.h5t_set_size(type_id, 1)
    API.h5t_set_cset(type_id, cset(C))
    Datatype(API.h5t_vlen_create(type_id))
end

# Opaque types
struct Opaque
  data
  tag::String
end

# An empty array type
struct EmptyArray{T} <: AbstractArray{T,0} end
# Required AbstractArray interface
Base.size(::EmptyArray) = ()
Base.IndexStyle(::Type{<:EmptyArray}) = IndexLinear()
Base.getindex(::EmptyArray, ::Int) = error("cannot index an `EmptyArray`")
Base.setindex!(::EmptyArray, v, ::Int) = error("cannot assign to an `EmptyArray`")
# Optional interface
Base.similar(::EmptyArray{T}) where {T} = EmptyArray{T}()
Base.similar(::EmptyArray, ::Type{S}) where {S} = EmptyArray{S}()
Base.similar(::EmptyArray, ::Type{S}, dims::Dims) where {S} = Array{S}(undef, dims)
# Override behavior for 0-dimensional Array
Base.length(::EmptyArray) = 0
# Required to avoid indexing during printing
Base.show(io::IO, E::EmptyArray) = print(io, typeof(E), "()")
Base.show(io::IO, ::MIME"text/plain", E::EmptyArray) = show(io, E)
# FIXME: Concatenation doesn't work for this type (it's treated as a length-1 array like
# Base's 0-dimensional arrays), so just forceably abort.
Base.cat_size(::EmptyArray) = error("concatenation of HDF5.EmptyArray is unsupported")
Base.cat_size(::EmptyArray, d) = error("concatenation of HDF5.EmptyArray is unsupported")

# Stub types to encode fixed-size arrays for API.H5T_ARRAY
struct FixedArray{T,D,L}
    data::NTuple{L,T}
end
Base.size(::Type{FixedArray{T,D,L}}) where {T,D,L} = D
Base.size(x::FixedArray) = size(typeof(x))
Base.eltype(::Type{FixedArray{T,D,L}}) where {T,D,L} = T
Base.eltype(x::FixedArray) = eltype(typeof(x))

struct FixedString{N,PAD}
    data::NTuple{N,UInt8}
end
Base.length(::Type{FixedString{N,PAD}}) where {N,PAD} = N
Base.length(str::FixedString) = length(typeof(str))
pad(::Type{FixedString{N,PAD}}) where {N,PAD} = PAD
pad(x::T) where {T<:FixedString} = pad(T)

struct VariableArray{T}
    len::Csize_t
    p::Ptr{Cvoid}
end
Base.eltype(::Type{VariableArray{T}}) where T = T

## Conversion between Julia types and HDF5 atomic types
hdf5_type_id(::Type{Bool})      = API.H5T_NATIVE_B8
hdf5_type_id(::Type{Int8})      = API.H5T_NATIVE_INT8
hdf5_type_id(::Type{UInt8})     = API.H5T_NATIVE_UINT8
hdf5_type_id(::Type{Int16})     = API.H5T_NATIVE_INT16
hdf5_type_id(::Type{UInt16})    = API.H5T_NATIVE_UINT16
hdf5_type_id(::Type{Int32})     = API.H5T_NATIVE_INT32
hdf5_type_id(::Type{UInt32})    = API.H5T_NATIVE_UINT32
hdf5_type_id(::Type{Int64})     = API.H5T_NATIVE_INT64
hdf5_type_id(::Type{UInt64})    = API.H5T_NATIVE_UINT64
hdf5_type_id(::Type{Float32})   = API.H5T_NATIVE_FLOAT
hdf5_type_id(::Type{Float64})   = API.H5T_NATIVE_DOUBLE
hdf5_type_id(::Type{Reference}) = API.H5T_STD_REF_OBJ
hdf5_type_id(::Type{<:AbstractString}) = API.H5T_C_S1

# It's not safe to use particular id codes because these can change, so we use characteristics of the type.
function _hdf5_type_map(class_id, is_signed, native_size)
    if class_id == API.H5T_INTEGER
        if is_signed == API.H5T_SGN_2
            return native_size == 1 ? Int8 :
                   native_size == 2 ? Int16 :
                   native_size == 4 ? Int32 :
                   native_size == 8 ? Int64 :
                   throw(KeyError((class_id, is_signed, native_size)))
        else
            return native_size == 1 ? UInt8 :
                   native_size == 2 ? UInt16 :
                   native_size == 4 ? UInt32 :
                   native_size == 8 ? UInt64 :
                   throw(KeyError((class_id, is_signed, native_size)))
        end
    else
        return native_size == 4 ? Float32 :
               native_size == 8 ? Float64 :
               throw(KeyError((class_id, is_signed, native_size)))
    end
end

# global configuration for complex support
const COMPLEX_SUPPORT = Ref(true)
const COMPLEX_FIELD_NAMES = Ref(("r", "i"))
enable_complex_support() = COMPLEX_SUPPORT[] = true
disable_complex_support() = COMPLEX_SUPPORT[] = false
set_complex_field_names(real::AbstractString, imag::AbstractString) =  COMPLEX_FIELD_NAMES[] = ((real, imag))

# Create a datatype from in-memory types
datatype(x::ScalarType) = Datatype(hdf5_type_id(typeof(x)), false)
datatype(::Type{T}) where {T<:ScalarType} = Datatype(hdf5_type_id(T), false)
datatype(A::AbstractArray{T}) where {T<:ScalarType} = Datatype(hdf5_type_id(T), false)
function datatype(::Type{Complex{T}}) where {T<:ScalarType}
  COMPLEX_SUPPORT[] || error("complex support disabled. call HDF5.enable_complex_support() to enable")
  dtype = API.h5t_create(API.H5T_COMPOUND, 2*sizeof(T))
  API.h5t_insert(dtype, COMPLEX_FIELD_NAMES[][1], 0, hdf5_type_id(T))
  API.h5t_insert(dtype, COMPLEX_FIELD_NAMES[][2], sizeof(T), hdf5_type_id(T))
  return Datatype(dtype)
end
datatype(x::Complex{<:ScalarType}) = datatype(typeof(x))
datatype(A::AbstractArray{Complex{T}}) where {T<:ScalarType} = datatype(eltype(A))

function datatype(str::AbstractString)
    type_id = API.h5t_copy(hdf5_type_id(typeof(str)))
    API.h5t_set_size(type_id, max(sizeof(str), 1))
    API.h5t_set_cset(type_id, cset(typeof(str)))
    Datatype(type_id)
end
function datatype(::Array{S}) where {S<:AbstractString}
    type_id = API.h5t_copy(hdf5_type_id(S))
    API.h5t_set_size(type_id, API.H5T_VARIABLE)
    API.h5t_set_cset(type_id, cset(S))
    Datatype(type_id)
end

# conversions to Julia types
function get_jl_type(obj_type::Datatype)
  class_id = API.h5t_get_class(obj_type)
  if class_id == API.H5T_OPAQUE
      return Opaque
  else
      return get_mem_compatible_jl_type(obj_type)
  end
end

function get_jl_type(obj)
  dtype = datatype(obj)
  try
      return get_jl_type(dtype)
  finally
      close(dtype)
  end
end

function get_mem_compatible_jl_type(obj_type::Datatype)
  class_id = API.h5t_get_class(obj_type)
  if class_id == API.H5T_STRING
      if API.h5t_is_variable_str(obj_type)
          return Cstring
      else
          N = sizeof(obj_type)
          PAD = API.h5t_get_strpad(obj_type)
          return FixedString{N,PAD}
      end
  elseif class_id == API.H5T_INTEGER || class_id == API.H5T_FLOAT
      native_type = API.h5t_get_native_type(obj_type)
      try
          native_size = API.h5t_get_size(native_type)
          if class_id == API.H5T_INTEGER
              is_signed = API.h5t_get_sign(native_type)
          else
              is_signed = nothing
          end
          return _hdf5_type_map(class_id, is_signed, native_size)
      finally
          API.h5t_close(native_type)
      end
  elseif class_id == API.H5T_BITFIELD
      return Bool
  elseif class_id == API.H5T_ENUM
      super_type = API.h5t_get_super(obj_type)
      try
          native_type = API.h5t_get_native_type(super_type)
          try
              native_size = API.h5t_get_size(native_type)
              is_signed = API.h5t_get_sign(native_type)
              return _hdf5_type_map(API.H5T_INTEGER, is_signed, native_size)
          finally
              API.h5t_close(native_type)
          end
      finally
          API.h5t_close(super_type)
      end
  elseif class_id == API.H5T_REFERENCE
      # TODO update to use version 1.12 reference functions/types
      return Reference
  elseif class_id == API.H5T_OPAQUE
      # TODO: opaque objects should get their own fixed-size data type; punning like
      #       this permits recursively reading (i.e. compound data type containing an
      #       opaque field). Requires figuring out what to do about the tag...
      len = Int(API.h5t_get_size(obj_type))
      return FixedArray{UInt8, (len,), len}
  elseif class_id == API.H5T_VLEN
      superid = API.h5t_get_super(obj_type)
      return VariableArray{get_mem_compatible_jl_type(Datatype(superid))}
  elseif class_id == API.H5T_COMPOUND
      N = API.h5t_get_nmembers(obj_type)

      membernames = ntuple(N) do i
          API.h5t_get_member_name(obj_type, i-1)
      end

      membertypes = ntuple(N) do i
          dtype = Datatype(API.h5t_get_member_type(obj_type, i-1))
          return get_mem_compatible_jl_type(dtype)
      end

      # check if should be interpreted as complex
      iscomplex = COMPLEX_SUPPORT[] &&
                  N == 2 &&
                  (membernames == COMPLEX_FIELD_NAMES[]) &&
                  (membertypes[1] == membertypes[2]) &&
                  (membertypes[1] <: ScalarType)

      if iscomplex
          return Complex{membertypes[1]}
      else
          return NamedTuple{Symbol.(membernames), Tuple{membertypes...}}
      end
  elseif class_id == API.H5T_ARRAY
      dims = API.h5t_get_array_dims(obj_type)
      nd = length(dims)
      eltyp = Datatype(API.h5t_get_super(obj_type))
      elT = get_mem_compatible_jl_type(eltyp)
      dimsizes = ntuple(i -> Int(dims[nd-i+1]), nd)  # reverse order
      return FixedArray{elT, dimsizes, prod(dimsizes)}
  end
  error("Class id ", class_id, " is not yet supported")
end

# convert special types to native julia types
function normalize_types(::Type{T}, buf::AbstractMatrix{UInt8}) where {T}
  # First dimension spans bytes of a single element of type T --- (recursively) normalize
  # each range of bytes to final type, returning vector of normalized data.
  return [_normalize_types(T, view(buf, :, ind)) for ind in axes(buf, 2)]
end

# high-level description which should always work --- here, the buffer contains the bytes
# for exactly 1 element of an object of type T, so reinterpret the `UInt8` vector as a
# length-1 array of type `T` and extract the (only) element.
function _typed_load(::Type{T}, buf::AbstractVector{UInt8}) where {T}
  return @inbounds reinterpret(T, buf)[1]
end
# fast-path for common concrete types with simple layout (which should be nearly all cases)
function _typed_load(::Type{T}, buf::V) where {T, V <: Union{Vector{UInt8}, Base.FastContiguousSubArray{UInt8,1}}}
  dest = Ref{T}()
  GC.@preserve dest buf Base._memcpy!(unsafe_convert(Ptr{Cvoid}, dest), pointer(buf), sizeof(T))
  return dest[]
  # TODO: The above can maybe be replaced with
  #   return GC.@preserve buf unsafe_load(convert(Ptr{t}, pointer(buf)))
  # dependent on data elements being properly aligned for all datatypes, on all
  # platforms.
end

_normalize_types(::Type{T}, buf::AbstractVector{UInt8}) where {T} = _typed_load(T, buf)
function _normalize_types(::Type{T}, buf::AbstractVector{UInt8}) where {K, T <: NamedTuple{K}}
  # Compound data types do not necessarily have members of uniform size, so instead of
  # dim-1 => bytes of single element and dim-2 => over elements, just loop over exact
  # byte ranges within the provided buffer vector.
  nv = ntuple(length(K)) do ii
      elT = fieldtype(T, ii)
      off = fieldoffset(T, ii) % Int
      sub = view(buf, off .+ (1:sizeof(elT)))
      return _normalize_types(elT, sub)
  end
  return NamedTuple{K}(nv)
end
function _normalize_types(::Type{V}, buf::AbstractVector{UInt8}) where {T, V <: VariableArray{T}}
  va = _typed_load(V, buf)
  pbuf = unsafe_wrap(Array, convert(Ptr{UInt8}, va.p), (sizeof(T), Int(va.len)))
  if do_normalize(T)
      # If `T` a non-trivial type, recursively normalize the vlen buffer.
      return normalize_types(T, pbuf)
  else
      # Otherwise if `T` is simple type, directly reinterpret the vlen buffer.
      # (copy since libhdf5 will reclaim `pbuf = va.p` in `API.h5d_vlen_reclaim`)
      return copy(vec(reinterpret(T, pbuf)))
  end
end
function _normalize_types(::Type{F}, buf::AbstractVector{UInt8}) where {T, F <: FixedArray{T}}
  if do_normalize(T)
      # If `T` a non-trivial type, recursively normalize the buffer after reshaping to
      # matrix with dim-1 => bytes of single element and dim-2 => over elements.
      return reshape(normalize_types(T, reshape(buf, sizeof(T), :)), size(F)...)
  else
      # Otherwise, if `T` is simple type, directly reinterpret the array and reshape to
      # final dimensions. The copy ensures (a) the returned array is independent of
      # [potentially much larger] read() buffer, and (b) that the returned data is an
      # Array and not ReshapedArray of ReinterpretArray of SubArray of ...
      return copy(reshape(reinterpret(T, buf), size(F)...))
  end
end
_normalize_types(::Type{Cstring}, buf::AbstractVector{UInt8}) = unsafe_string(_typed_load(Ptr{UInt8}, buf))
_normalize_types(::Type{T}, buf::AbstractVector{UInt8}) where {T <: FixedString} = unpad(String(buf), pad(T))

do_normalize(::Type{T}) where {T} = false
do_normalize(::Type{NamedTuple{T,U}}) where {U,T} = any(i -> do_normalize(fieldtype(U,i)), 1:fieldcount(U))
do_normalize(::Type{T}) where {T <: Union{Cstring,FixedString,FixedArray,VariableArray}} = true

do_reclaim(::Type{T}) where {T} = false
do_reclaim(::Type{NamedTuple{T,U}}) where {U,T} = any(i -> do_reclaim(fieldtype(U,i)), 1:fieldcount(U))
do_reclaim(::Type{T}) where T <: Union{Cstring,VariableArray} = true


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

function dottest(n, iter)
    a = rand(n)
    b = rand(n)
    c = similar(a)
    for i = 1:iter
        c = dot(a,b)
    end
    c
end

function axpytest(n, iter)
    a = rand(1)[1]
    x = rand(n)
    y = zeros(n)
    for i = 1:iter
        BLAS.axpy!(a, x, y)
    end
    y
end

problemsizes = [(2, 10^6, "tiny"), (2^4, 10^6, "small"), (2^6, 10^6, "medium"), (2^8, 10^5, "large"), (2^10, 10^5, "huge")]
testdata = [(dottest, "dot", "dot product", problemsizes), (axpytest, "axpy", "axpy", problemsizes)]
include("../perfgeneric.jl")

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia