text:
function part1(input)
    seats = reduce(hcat, (line -> split(line, "") .== "L").(readlines(input)))
    occupied = falses(size(seats))
    while true
        neighbors = count_neighbors(occupied)
        next_occupied = ((neighbors .== 0) .|
                         ((neighbors .< 4) .& occupied)) .&
                         seats
        occupied == next_occupied && return count(occupied)
        occupied = next_occupied
    end
end

function part2(input)
    seats = reduce(hcat, (line -> split(line, "") .== "L").(readlines(input)))
    occupied = falses(size(seats))
    while true
        neighbors = count_visible(occupied, seats)
        next_occupied = ((neighbors .== 0) .|
                         ((neighbors .< 5) .& occupied)) .&
                         seats
        occupied == next_occupied && return count(occupied)
        occupied = next_occupied
    end
end

function count_neighbors(occupied)
    neighbors = zeros(Int, size(occupied))
    indices = CartesianIndices(occupied)
    for index in indices
        for Δ in CartesianIndices((-1:1, -1:1))
            Δ == CartesianIndex(0, 0) && continue
            if index + Δ in indices
                neighbors[index] += occupied[index + Δ]
            end
        end
    end
    return neighbors
end

function count_visible(occupied, seats)
    visible = zeros(Int, size(occupied))
    indices = CartesianIndices(occupied)
    for index in indices
        seats[index] || continue
        for Δ in CartesianIndices((-1:1, -1:1))
            Δ == CartesianIndex(0, 0) && continue
            current_index = index + Δ
            while current_index in indices
                if seats[current_index]
                    visible[index] += occupied[current_index]
                    break
                end
                current_index += Δ
            end
        end
    end
    return visible
end


text:
"""
    add(x, y)

Perform addition between two floating point numbers.

# Arguments
- `a::Float64`: Number 1.
- `b::Float64`: Number 2.

# Returns
- `a + b::Float64`

# Examples
```julia-repl
julia> add(1.0, 2.0)
3.0
```
"""
add(a::Float64, b::Float64) = a + b

ext: jl


text:
#Tests various bifurcation capabilities.
for rn in reaction_networks_standard
    p_vals = map(amp->amp*rand(length(rn.params)),[1.,10.,100.])
    param1 = rand(rn.params)
    param2 = rand(rn.params)
    range_tupple = (1.,rand(5.:10.))
    range1 = 1.: rand(5.:10.)
    range2 = 1.: rand(5.:10.)
    for p in p_vals
        bif1 = bifurcations(rn,p,param1,range_tupple)
        bif2 = bifurcations(rn,SimpleHCBifurcationSolver(),p,param1,range_tupple)
        bif_grid = bifurcation_grid(rn,p,param1,range1)
        bif_grid_2d = bifurcation_grid_2d(rn,p,param1,range1,param2,range2)
        bif_grid_dia = bifurcation_grid_diagram(rn,p,param1,range1,param2,range_tupple)
        plot(bif1); plot(bif2); plot(bif_grid); plot(bif_grid_2d); plot(bif_grid_dia);
    end
end

#More extensive tests on familiar networks.
brusselator_network = @reaction_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;
brusselator_p = [1.,4.]
@make_hc_template brusselator_network
@add_hc_template brusselator_network
bif_brusselator = bifurcations(brusselator_network,brusselator_p,:B,(1.,4.))
if VERSION >=  v"1.1"
    @test length(bif_brusselator.paths) == 1
else
    @warn "Some test for bifurcation diagram path numbers are disabled for Julia 1.0 versions. Bifurcation diagrams might contain slight errors for tehse versions."
end

σ_network = @reaction_network begin
    v0 + hill(σ,v,K,n), ∅ → (σ+A)
    deg, (σ,A,Aσ) → ∅
    (kB,kD), A + σ ↔ Aσ
    S*kC, Aσ → σ
end v0 v K n kD kB kC deg S;
fix_parameters(σ_network,n=4)
σ_p = [0.005, 0.1, 2.8, 4, 10, 100, 0.1, 0.01, 0.5]
make_hc_template!(σ_network)
add_hc_template!(σ_network)
ss_σ = steady_states(σ_network,σ_p)
@test(length(ss_σ)==3)
stabs_σ = stability(ss_σ,σ_network,σ_p)
@test(sum(stabs_σ.==true)==2)
bif_σ = bifurcations(σ_network,σ_p,:S,(0.,2.))
if VERSION >= v"1.1"
    @test length(bif_σ.paths) == 3
else
    @warn "Some test for bifurcation diagram path numbers are disabled for Julia 1.0 versions. Bifurcation diagrams might contain slight errors for tehse versions."
end

cc_network = @reaction_network begin
  k1, 0 --> Y
  k2p, Y --> 0
  k2pp*P, Y --> 0
  (k3p+k3pp*A)/(J3+Po), Po-->P
  (k4*m)/(J4+P), Y + P --> Y + Po
end k1 k2p k2pp k3p k3pp A J3 k4 m J4
cc_p = [0.04,0.04,1.,1.,10.0,0.,0.04,35.,.3,.04]
@add_constraint cc_network P+Po=1
add_hc_template!(cc_network)
@add_hc_template cc_network
ss_cc = steady_states(cc_network,cc_p)
@test(length(ss_cc)==3)
stabs_cc = stability(ss_σ,σ_network,σ_p)
@test(sum(stabs_cc.==true)==2)
bif_cc = bifurcations(cc_network,cc_p,:m,(.01,.65))
if VERSION>= v"1.1"
    @test length(bif_cc.paths) == 3
else
    @warn "Some test for bifurcation diagram path numbers are disabled for Julia 1.0 versions. Bifurcation diagrams might contain slight errors for tehse versions."
end

bs_network = @reaction_network begin
    d,    (X,Y) → ∅
    hillR(Y,v1,K1,n1), ∅ → X
    hillR(X,v2,K2,n2), ∅ → Y
end d v1 K1 n1 v2 K2 n2
bs_p = [0.01, 1. , 30., 3, 1., 30, 3];
bif_bs = bifurcations(bs_network, bs_p,:v1,(.1,10.))
if VERSION >= v"1.1"
    @test length(bif_bs.paths) == 3
    @test sum(map(p->median(p.stability_types),bif_bs.paths))==2
else
    @warn "Some test for bifurcation diagram path numbers are disabled for Julia 1.0 versions. Bifurcation diagrams might contain slight errors for tehse versions."
end


text:
include("Parsers.jl")
include("Utils.jl")
using DataInterpolations

abstract type AbstractDiff  end

struct AnalyticalDeriv <: AbstractDiff end
struct FiniteDiff <: AbstractDiff end
struct TotalVariationalDerivativative <: AbstractDiff end

function differentiate(sol,o::TotalVariationalDerivativative)
    # interp = LagrangeInterpolation(munge(sol.u),sol.t)
    interp = QuadraticInterpolation(munge(sol.u),sol.t)
    y_reconstruct = interp.(sol.t)
    v = [DataInterpolations.derivative(interp,k) for k in sol.t]
    w = munge(v)
    if length(sol.u[1]) == 3
        fig = plot(w[1,:],w[2,:],w[3,:])
    elseif length(sol.u[1]) == 2
        fig = plot(w[1,:],w[2,:])
    end
    display(fig)
    return v
end

function differentiate(sol,o::FiniteDiff)
    q = sol.u[2:end] .- sol.u[1:end-1]
    q = q./(sol.t[2] - sol.t[1])
    w = munge(q)
    if length(sol.u[1]) == 3
        fig = plot(w[1,:],w[2,:],w[3,:])
    elseif length(sol.u[1]) == 2
        fig = plot(w[1,:],w[2,:])
    end
    display(fig)
    return w
end

function differentiate(sol::Matrix,t,o::TotalVariationalDerivativative)
    interp = QuadraticInterpolation(sol,t)
    y_reconstruct = interp.(t)
    v = [DataInterpolations.derivative(interp,k) for k in t]
    w = munge(v)
    if size(sol,1) == 3
        fig = plot(w[1,:],w[2,:],w[3,:])
    elseif size(sol,1) == 2
        fig = plot(w[1,:],w[2,:])
    end
    display(fig)
    return v
end

function differentiate(sol::Matrix,t,o::FiniteDiff)
    q = sol[:,2:end] .- sol[:,1:end-1]
    q = q./(t[2] - t[1])
    w = q
    if size(sol,1) == 3
        fig = plot(w[1,:],w[2,:],w[3,:])
    elseif size(sol,1) == 2
        fig = plot(w[1,:],w[2,:])
    end
    display(fig)
    return w
end

#--------------- The analytical derivative is computed for testing---------------------#

function differentiate(sol,ds::LorenzSystem,o::AnalyticalDeriv)
    function lorenz(u) # Modify when changing parameters
        u1 = 10*(u[2] - u[1])
        u2 = u[1]*(28 - u[3]) - u[2]
        u3 = u[1]*u[2] - 8/3 * u[3]
        return [u1 u2 u3]
    end
    deriv = []
    for i = 1:length(sol.u)
        append!(deriv,lorenz(sol.u[i]))
    end
    deriv = reshape(deriv,(3,length(sol.u)))
    fig = plot(deriv[1,:],deriv[2,:],deriv[3,:],title="Lorenz - Derivative")
    display(fig)
    return deriv

end

function differentiate(sol,ds::LotkaVolterra,o::AnalyticalDeriv)
    function lk(u) # Modify when changing parameters
        u1 = 0.7*u[1] - 0.3*u[1]*u[2]
        u2 = -0.3*u[2] + 0.4*u[1]*u[2]
        return [u1 u2]
    end
    deriv = []
    for i = 1:length(sol.u)
        append!(deriv,lk(sol.u[i]))
    end
    deriv = reshape(deriv,(2,length(sol.u)))
    fig = plot(deriv[1,:],deriv[2,:],title="Lotka Volterra - Derivative")
    display(fig)
    return deriv
end


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

module Pkg2Types

export VersionInterval, VersionSet
import Base: show, isempty, in, intersect, union!, union, ==, hash, copy, deepcopy_internal

import ...Types: VersionBound, VersionRange

struct VersionInterval
    lower::VersionNumber
    upper::VersionNumber
end
VersionInterval(lower::VersionNumber) = VersionInterval(lower,typemax(VersionNumber))
VersionInterval() = VersionInterval(typemin(VersionNumber))

show(io::IO, i::VersionInterval) = print(io, "[$(i.lower),$(i.upper))")
isempty(i::VersionInterval) = i.upper <= i.lower
in(v::VersionNumber, i::VersionInterval) = i.lower <= v < i.upper
intersect(a::VersionInterval, b::VersionInterval) = VersionInterval(max(a.lower,b.lower), min(a.upper,b.upper))
==(a::VersionInterval, b::VersionInterval) = a.lower == b.lower && a.upper == b.upper
hash(i::VersionInterval, h::UInt) = hash((i.lower, i.upper), h + (0x0f870a92db508386 % UInt))

# Used to translate from Pkg2 intervals to Pkg ranges
function Base.convert(::Type{VersionRange}, a::VersionInterval)
    lower, upper = a.lower, a.upper

    lb = VersionBound(lower.major, lower.minor, lower.patch)

    vb = UInt32[upper.major, upper.minor, upper.patch]
    # Consider e.g. "0.1.0 0.1.0+" to be the same as VerionRange("0.1.0")
    # by making the "+" bump the patch version
    if typeof(upper.build) == Tuple{String} && vb[3] != typemax(UInt32)
        vb[3] += 1
    end
    i = 3
    while i > 0 && vb[i] == 0
        pop!(vb)
        i -= 1
    end
    # NOTE: an upper bound of 0 could happen in principle, e.g. [v"0.0.0-", v"0.0.0")
    #       but we just ignore this here...
    i > 0 || error("invalid interval upper bound v$upper")
    vb[i] -= 1
    ub = VersionBound(vb...)

    return VersionRange(lb, ub)
end

function normalize!(ivals::Vector{VersionInterval})
    # VersionSet internal normalization:
    # removes empty intervals and fuses intervals without gaps
    # e.g.:
    #     [0.0.0,1.0.0) ∪ [1.0.0,1.5.0) ∪ [1.6.0,1.6.0) ∪ [2.0.0,∞)
    # becomes:
    #     [0.0.0,1.5.0) ∪ [2.0.0,∞)
    # (still assumes that lower bounds are sorted, and intervals do
    # not overlap)
    l = length(ivals)
    l == 0 && return ivals

    lo, up, k0 = ivals[1].lower, ivals[1].upper, 1
    fusing = false
    for k = 2:l
        lo1, up1 = ivals[k].lower, ivals[k].upper
        if lo1 == up
            up = up1
            fusing = true
            continue
        end
        if lo < up
            # The only purpose of the "fusing" check is to avoid
            # extra allocations
            ivals[k0] = fusing ? VersionInterval(lo, up) : ivals[k-1]
            k0 += 1
        end
        fusing = false
        lo, up = lo1, up1
    end
    if lo < up
        ivals[k0] = fusing ? VersionInterval(lo, up) : ivals[l]
        k0 += 1
    end
    resize!(ivals, k0 - 1)
    return ivals
end

struct VersionSet
    intervals::Vector{VersionInterval}
    VersionSet(intervals::Vector{VersionInterval}) = new(normalize!(intervals))
    # copy is defined inside the struct block to call `new` directly
    # without going through `normalize!`
    Base.copy(vset::VersionSet) = new(copy(vset.intervals))
end
function VersionSet(versions::Vector{VersionNumber})
    intervals = VersionInterval[]
    if isempty(versions)
        push!(intervals, VersionInterval())
    else
        isodd(length(versions)) && push!(versions, typemax(VersionNumber))
        while !isempty(versions)
            push!(intervals, VersionInterval(popfirst!(versions), popfirst!(versions)))
        end
    end
    VersionSet(intervals)
end
VersionSet(versions::VersionNumber...) = VersionSet(VersionNumber[versions...])

const empty_versionset = VersionSet(VersionInterval[])

# Windows console doesn't like Unicode
const _empty_symbol = @static Sys.iswindows() ? "empty" : "∅"
const _union_symbol = @static Sys.iswindows() ? " or " : " ∪ "
show(io::IO, s::VersionSet) = isempty(s) ? print(io, _empty_symbol) :
                                           join(io, s.intervals, _union_symbol)
isempty(s::VersionSet) = all(isempty, s.intervals)
in(v::VersionNumber, s::VersionSet) = any(i->in(v,i), s.intervals)
function intersect(A::VersionSet, B::VersionSet)
    (isempty(A) || isempty(B)) && return copy(empty_versionset)
    ivals = [intersect(a,b) for a in A.intervals for b in B.intervals]
    sort!(ivals, by=i->i.lower)
    VersionSet(ivals)
end

union(A::VersionSet, B::VersionSet) = union!(copy(A), B)
function union!(A::VersionSet, B::VersionSet)
    A == B && return A
    ivals = A.intervals
    for intB in B.intervals
        lB, uB = intB.lower, intB.upper
        k0 = findfirst(i->(i.upper > lB), ivals)
        if k0 == 0
            push!(ivals, intB)
            continue
        end
        lB = min(lB, ivals[k0].lower)
        for k1 = k0:length(ivals)
            intA = ivals[k1]
            if uB < intA.lower
                splice!(ivals, k0:(k1-1), (VersionInterval(lB, uB),))
                break
            elseif uB ∈ intA || k1 == length(ivals)
                splice!(ivals, k0:k1, (VersionInterval(lB, max(uB, intA.upper)),))
                break
            end
        end
    end
    normalize!(ivals)
    return A
end

==(A::VersionSet, B::VersionSet) = A.intervals == B.intervals
hash(s::VersionSet, h::UInt) = hash(s.intervals, h + (0x2fd2ca6efa023f44 % UInt))
deepcopy_internal(vs::VersionSet, ::IdDict) = copy(vs)

end # module

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
function lower_load_scalar!(
    q::Expr, op::Operation, ua::UnrollArgs, umin::Int, inds_calc_by_ptr_offset::Vector{Bool}
)
    loopdeps = loopdependencies(op)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    @assert vectorized ∉ loopdeps
    # mvar, opu₁, opu₂ = variable_name_and_unrolled(op, u₁loop, u₂loop, suffix)
    mvar = variable_name(op, suffix)
    opu₁ = isu₁unrolled(op)
    ptr = vptr(op)
    U = opu₁ ? u₁ : 1
    if instruction(op).instr !== :conditionalload
        for u ∈ umin:U-1
            varname = varassignname(mvar, u, opu₁)
            td = UnrollArgs(ua, u)
            push!(q.args, Expr(:(=), varname, Expr(:call, lv(:vload), ptr, mem_offset_u(op, td, inds_calc_by_ptr_offset))))
        end
    else
        opu₂ = !isnothing(suffix) && u₂loopsym ∈ loopdeps
        condop = last(parents(op))
        condvar, condu₁ = condvarname_and_unroll(condop, u₁loopsym, u₂loopsym, suffix, opu₂)
        for u ∈ umin:U-1
            condsym = varassignname(condvar, u, condu₁)
            varname = varassignname(mvar, u, u₁loopsym ∈ loopdependencies(op))
            td = UnrollArgs(ua, u)
            load = Expr(:call, lv(:vload), ptr, mem_offset_u(op, td, inds_calc_by_ptr_offset))
            cload = Expr(:if, condsym, load, Expr(:call, :zero, Expr(:call, :eltype, ptr)))
            push!(q.args, Expr(:(=), varname, cload))
        end
    end
    nothing
end
function pushvectorload!(
    q::Expr, op::Operation, var::Symbol, td::UnrollArgs, U::Int, vectorized::Symbol, mask, u₁unrolled::Bool, inds_calc_by_ptr_offset::Vector{Bool}
)
    @unpack u₁, u₁loopsym, u₂loopsym, suffix = td
    ptr = vptr(op)
    vecnotunrolled = vectorized !== u₁loopsym
    name, mo = name_memoffset(var, op, td, u₁unrolled, inds_calc_by_ptr_offset)
    instrcall = Expr(:call, lv(:vload), ptr, mo)

    iscondstore = instruction(op).instr === :conditionalload
    maskend = mask !== nothing && (vecnotunrolled || u₁ == U - 1)
    if iscondstore
        condop = last(parents(op))
        opu₂ = !isnothing(suffix) && u₂loopsym ∈ loopdependencies(op)
        condvar, condu₁ = condvarname_and_unroll(condop, u₁loopsym, u₂loopsym, suffix, opu₂)
        condsym = varassignname(condvar, u₁, condu₁)
        if vectorized ∈ loopdependencies(condop)
            if maskend
                push!(instrcall.args, Expr(:call, :&, condsym, mask))
            else
                push!(instrcall.args, condsym)
            end
        else
            if maskend
                push!(instrcall.args, mask)
            end
            instrcall = Expr(:if, condsym, instrcall, Expr(:call, lv(:vzero), VECTORWIDTHSYMBOL, Expr(:call, :eltype, ptr)))
        end
    elseif maskend
        push!(instrcall.args, mask)
    end
    push!(q.args, Expr(:(=), name, instrcall))
    # push!(q.args, :(@show $name))
end
function prefetchisagoodidea(ls::LoopSet, op::Operation, td::UnrollArgs)
    # return false
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, u₂max, suffix = td
    vectorized ∈ loopdependencies(op) || return 0
    u₂loopsym === Symbol("##undefined##") && return 0
    dontskip = (VectorizationBase.CACHELINE_SIZE ÷ VectorizationBase.REGISTER_SIZE) - 1
    # u₂loopsym is vectorized
    # u₁vectorized = vectorized === u₁loopsym
    u₂vectorized = vectorized === u₂loopsym
    (!isnothing(suffix) && u₂vectorized && !iszero(suffix & dontskip)) && return 0
    innermostloopsym = first(names(ls))
    loopedindex = op.ref.loopedindex
    if length(loopedindex) > 1 && first(loopedindex)
        indices = getindices(op)
        if first(indices) === vectorized && last(indices) === innermostloopsym
            # We want at least 4 reuses per load
            uses = ifelse(isu₁unrolled(op), 1, u₁)
            uses = ifelse(isu₂unrolled(op), uses, uses * u₂max)
            uses < 4 && return 0
            innermostloopindv = findall(map(isequal(innermostloopsym), getindices(op)))
            isone(length(innermostloopindv)) || return 0
            innermostloopind = first(innermostloopindv)
            if prod(s -> Float64(length(getloop(ls, s))), @view(indices[1:innermostloopind-1])) ≥ 120.0 && length(getloop(ls, innermostloopsym)) ≥ 120
                if op.ref.ref.offsets[innermostloopind] < 120
                    for opp ∈ operations(ls)
                        iscompute(opp) && (innermostloopsym ∈ loopdependencies(opp)) && load_constrained(opp, u₁loopsym, u₂loopsym, true) && return 0
                    end
                    return innermostloopind
                end
            end
        end
    end
    0
end
function add_prefetches!(q::Expr, ls::LoopSet, op::Operation, td::UnrollArgs, prefetchind::Int, umin::Int)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, u₂max = td
    dontskip = (64 ÷ VectorizationBase.REGISTER_SIZE) - 1
    ptr = vptr(op)
    innermostloopsym = first(names(ls))
    us = ls.unrollspecification[]
    prefetch_distance = u₁loopsym === innermostloopsym ? us.u₁ : ( u₂loopsym === innermostloopsym ? us.u₂ : 1 )
    # prefetch_distance = u₁loopsym === innermostloopsym ? u₁ : ( u₂loopsym === innermostloopsym ? u₂max : 1 )
    prefetch_multiplier = 5
    prefetch_distance *= prefetch_multiplier
    offsets = op.ref.ref.offsets
    inner_offset = offsets[prefetchind]
    offsets[prefetchind] = inner_offset + prefetch_distance
    gespinds = mem_offset_u(op, UnrollArgs(td, 0), indices_calculated_by_pointer_offsets(ls, op.ref))
    offsets[prefetchind] = inner_offset
    ptr = vptr(op)
    gptr = Symbol(ptr, "##GESPEDPREFETCH##")
    for i ∈ eachindex(gespinds.args)
       gespinds.args[i] = Expr(:call, lv(:extract_data), gespinds.args[i])
    end
    push!(q.args, Expr(:(=), gptr, Expr(:call, lv(:gesp), ptr, gespinds)))

    inds = Expr(:tuple)
    indices = getindicesonly(op)

    i = 0
    for (j,ind) ∈ enumerate(indices)
        push!(inds.args, Expr(:call, lv(:Zero)))
        (ind == u₁loopsym) && (i = j)
    end
    push!(q.args, Expr(:call, lv(:prefetch0), gptr, copy(inds)))

    for u ∈ 1+umin:u₁
        # for u ∈ umin:min(umin,U-1)
        (u₁loopsym === vectorized && !iszero(u & dontskip)) && continue
        if u₁loopsym === vectorized
            if isone(u)
                inds.args[i] = Expr(:call, lv(:unwrap), VECTORWIDTHSYMBOL)
            else
                inds.args[i] = Expr(:call, lv(:valmul), VECTORWIDTHSYMBOL, u)
            end
        else
            inds.args[i] = Expr(:call, Expr(:curly, lv(:Static), u))
        end
        push!(q.args, Expr(:call, lv(:prefetch0), gptr, copy(inds)))
    end
    nothing
end
function lower_load_vectorized!(
    q::Expr, ls::LoopSet, op::Operation, td::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned} = nothing, umin::Int = 0
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = td
    loopdeps = loopdependencies(op)
    @assert isvectorized(op)
    opu₁ = isu₁unrolled(op)
    inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
    if opu₁
        umin = umin
        U = u₁
    else
        umin = -1
        U = 0
    end
    # Urange = unrolled ∈ loopdeps ? 0:U-1 : 0
    var = variable_name(op, suffix)
    for u ∈ umin:U-1
        td = UnrollArgs(td, u)
        pushvectorload!(q, op, var, td, U, vectorized, mask, opu₁, inds_calc_by_ptr_offset)
    end
    prefetchind = prefetchisagoodidea(ls, op, td)
    iszero(prefetchind) || add_prefetches!(q, ls, op, td, prefetchind, umin)
    nothing
end
function indisvectorized(ls::LoopSet, ind::Symbol, vectorized::Symbol)
    for op ∈ operations(ls)
        ((op.variable === ind) && isvectorized(op)) && return true
    end
    false
end

# function lower_load_for_optranslation!(
#     q::Expr, op::Operation, ls::LoopSet, td::UnrollArgs{Int}, mask::Union{Nothing,Symbol,Unsigned}, translationind::Int
# )
#     @unpack u₁, u₁loopsym, u₂loopsym, vectorized, u₂max, suffix = td
#     iszero(suffix) || return

#     gespinds = mem_offset(op, UnrollArgs(td, 0), indices_calculated_by_pointer_offsets(ls, op.ref), false)
#     ptr = vptr(op)
#     gptr = Symbol(ptr, "##GESPED##")
#     for i ∈ eachindex(gespinds.args)
#         if i != translationind
#             gespinds.args[i] = Expr(:call, lv(:extract_data), gespinds.args[i])
#         end
#     end
#     push!(q.args, Expr(:(=), gptr, Expr(:call, lv(:gesp), ptr, gespinds)))

#     inds = Expr(:tuple)
#     ginds = Expr(:tuple)
#     indices = getindicesonly(op)

#     for (i,ind) ∈ enumerate(indices)
#         if i == translationind # ind cannot be the translation ind
#             push!(inds.args, Expr(:call, lv(:Zero)))
#             push!(ginds.args, Expr(:call, Expr(:curly, lv(:Static), 1)))
#         elseif (ind === vectorized) || indisvectorized(ls, ind, vectorized)
#             push!(inds.args, _MMind(Expr(:call, lv(:Zero))))
#             push!(ginds.args, Expr(:call, lv(:Zero)))
#         else
#             push!(inds.args, Expr(:call, lv(:Zero)))
#             push!(ginds.args, Expr(:call, lv(:Zero)))
#         end
#     end
#     varbase = variable_name(op, 0)
#     vloadexpr = Expr(:call, lv(:vload), gptr, inds)
#     gespexpr  = Expr(:(=), gptr, Expr(:call, lv(:gesp), gptr, ginds))
#     push!(q.args, Expr(:(=), Symbol(varbase, 0), vloadexpr))

#     for u ∈ 1:u₁-1
#         push!(q.args, gespexpr)
#         push!(q.args, Expr(:(=), Symbol(varbase, u), vloadexpr))
#     end
#     # this takes care of u₂ == 0
#     offset = u₁
#     for u₂ ∈ 1:u₂max-1
#         varold = varbase
#         varbase = variable_name(op, u₂)
#         for u ∈ 0:u₁-2
#             push!(q.args, Expr(:(=), Symbol(varbase, u), Symbol(varold, u + 1)))
#         end
#         push!(q.args, gespexpr)
#         push!(q.args, Expr(:(=), Symbol(varbase, u₁ - 1), vloadexpr))
#         offset += 1
#     end
#     nothing
# end

function lower_load_for_optranslation!(
    q::Expr, op::Operation, ls::LoopSet, td::UnrollArgs{Int}, mask::Union{Nothing,Symbol,Unsigned}, translationind::Int
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, u₂max, suffix = td
    iszero(suffix) || return

    gespinds = mem_offset(op, UnrollArgs(td, 0), indices_calculated_by_pointer_offsets(ls, op.ref), false)
    ptr = vptr(op)
    gptr = Symbol(ptr, "##GESPED##")
    for i ∈ eachindex(gespinds.args)
        if i != translationind
            gespinds.args[i] = Expr(:call, lv(:extract_data), gespinds.args[i])
        end
    end
    push!(q.args, Expr(:(=), gptr, Expr(:call, lv(:gesp), ptr, gespinds)))

    inds = Expr(:tuple)
    indices = getindicesonly(op)

    for (i,ind) ∈ enumerate(indices)
        if i == translationind # ind cannot be the translation ind
            push!(inds.args, Expr(:call, Expr(:curly, lv(:Static), 0)))
        elseif (ind === vectorized) || indisvectorized(ls, ind, vectorized)
            push!(inds.args, _MMind(Expr(:call, lv(:Zero))))
        else
            push!(inds.args, Expr(:call, lv(:Zero)))
        end
    end
    varbase = variable_name(op, 0)
    push!(q.args, Expr(:(=), Symbol(varbase, 0), Expr(:call, lv(:vload), gptr, copy(inds))))

    for u ∈ 1:u₁-1
        inds.args[translationind] = Expr(:call, Expr(:curly, lv(:Static), u))
        push!(q.args, Expr(:(=), Symbol(varbase, u), Expr(:call, lv(:vload), gptr, copy(inds))))
    end
    # this takes care of u₂ == 0
    offset = u₁
    for u₂ ∈ 1:u₂max-1
        varold = varbase
        varbase = variable_name(op, u₂)
        for u ∈ 0:u₁-2
            push!(q.args, Expr(:(=), Symbol(varbase, u), Symbol(varold, u + 1)))
        end
        inds.args[translationind] = Expr(:call, Expr(:curly, lv(:Static), offset))
        push!(q.args, Expr(:(=), Symbol(varbase, u₁ - 1), Expr(:call, lv(:vload), gptr, copy(inds))))
        offset += 1
    end
    nothing
end

# TODO: this code should be rewritten to be more "orthogonal", so that we're just combining separate pieces.
# Using sentinel values (eg, T = -1 for non tiling) in part to avoid recompilation.
function lower_load!(
    q::Expr, op::Operation, ls::LoopSet, td::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned} = nothing
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = td
    if !isnothing(suffix) && ls.loadelimination[]
        istr, ispl = isoptranslation(ls, op, UnrollSymbols(u₁loopsym, u₂loopsym, vectorized))
        if !iszero(istr) & ispl
            return lower_load_for_optranslation!(q, op, ls, td, mask, istr)
        elseif suffix > 0
            if u₂loopsym !== vectorized
                mno, id = maxnegativeoffset(ls, op, u₂loopsym)
                if -suffix < mno < 0
                    varnew = variable_name(op, suffix)
                    varold = variable_name(operations(ls)[id], suffix + mno)
                    opold = operations(ls)[id]
                    if isu₁unrolled(op)
                        for u ∈ 0:u₁-1
                            push!(q.args, Expr(:(=), Symbol(varnew, u), Symbol(varold, u)))
                        end
                    else
                        push!(q.args, Expr(:(=), varnew, varold))
                    end
                    return
                end
            end
        end
    end
    umin = 0
    if isvectorized(op)
        lower_load_vectorized!(q, ls, op, td, mask, umin)
    else
        lower_load_scalar!(q, op, td, umin, indices_calculated_by_pointer_offsets(ls, op.ref))
    end
end


text:
abstract type AbstractPlant <: AbstractEnv end

struct DelayFDI_Plant <: AbstractPlant
    multicopter::Multicopter
    fdi::DelayFDI
    faults::Vector{AbstractFault}
end

function State(env::DelayFDI_Plant)
    @unpack multicopter, fdi, faults = env
    return function (; args_multicopter=())
        x_multicopter = State(multicopter)(args_multicopter...)
        ComponentArray(multicopter=x_multicopter)
    end
end

function EffectivenessMatrixFunction(env::DelayFDI_Plant)
    @unpack multicopter, fdi, faults = env
    @unpack dim_input = multicopter
    # actuator faults
    actuator_faults = faults |> Filter(fault -> typeof(fault) <: AbstractActuatorFault) |> collect
    affect_actuator! = Affect!(AbstractActuatorFault[actuator_faults...])
    # effectiveness matrix
    Λ_func = function (t)
        Λ = ones(dim_input) |> Diagonal
        affect_actuator!(Λ, t)
        Λ
    end
    # delayed effectiveness matrix estimation
    Λ̂_func = Estimate(fdi, Λ_func)
    (; Λ_func=Λ_func, Λ̂_func=Λ̂_func)
end

function Dynamics!(env::DelayFDI_Plant)
    @unpack multicopter, fdi, faults = env
    effectiveness_matrix_functions = EffectivenessMatrixFunction(env)
    @unpack Λ_func, Λ̂_func = effectiveness_matrix_functions
    @Loggable function dynamics!(dX, X, p, t; u)
        Λ = Λ_func(t)
        @nested_log :FDI Λ̂ = Λ̂_func(t)
        # @nested_log :input u_saturated = pinv(Λ̂) * FlightSims.saturate(multicopter, Λ̂*u)  # saturation considering fault
        # @nested_log :input u_saturated = u  # saturation considering fault
        # @nested_log Dynamics!(multicopter)(dX.multicopter, X.multicopter, (), t; u=u_saturated, Λ=Λ)
        @nested_log Dynamics!(multicopter)(dX.multicopter, X.multicopter, (), t; u=u, Λ=Λ)
    end
end


text:
using SurfaceGeometry
using LinearAlgebra
using Test

# Pkg.test("SurfaceGeometry")
# write your own tests here
# @test 1 == 2

# Loading of spherical mesh

function SphereError(points)
    error = 0
    for xkey in 1:size(points,2)
        pos = points[:,xkey]
        n = pos/norm(pos)
        err = norm(n - pos)
        if err>error
            error = err
        end
    end
    return error
end

### Loading of mesh for other tests
@info "Mesh loading test"
using JLD
data = load(joinpath(dirname(@__FILE__),"sphere.jld"))
points = data["points"]
faces = data["faces"]

@test SphereError(points) < 0.01

@info "Topology function tests"

t = [5,7,10,7,5,6,4,0,3,0,4,6,4,7,6,4,9,10,7,4,10,0,2,1,2,0,6,2,5,1,5,2,6,8,4,3,4,11,9,8,11,4,9,11,3,11,8,3]
faces = reshape(t,(3,div(length(t),3))) .+ 1

triangles = []
for i in FaceVRing(5,faces)
    #println("i")
    push!(triangles,i)
end

@test sort(triangles)==[3,4,5,6,7,12,13,14]

triverticies = []
for i in DoubleVertexVRing(5,faces)
    #println("i")
    push!(triverticies,i)
end

@test (1,4) in triverticies

verticies = []
for i in VertexVRing(5,faces)
    push!(verticies,i)
end

@test sort(verticies)==[1,4,7,8,9,10,11,12] #[1,2,6,7]

@info "Testing Complex DS"

points = zero(faces)
fb = FaceBasedDS(faces)

triangles = []
for i in FaceVRing(5,fb)
    push!(triangles,i)
end
@test sort(triangles)==[3,4,5,6,7,12,13,14]

triverticies = []
for i in DoubleVertexVRing(5,fb)
    push!(triverticies,i)
end

@test (1,4) in triverticies

verticies = []
for i in VertexVRing(5,fb)
     push!(verticies,i)
end

@test sort(verticies)==[1,4,7,8,9,10,11,12]

@info "Topology tests for Connectivity table"

using JLD
data = load(joinpath(dirname(@__FILE__),"sphere.jld"))
#data = load("sphere.jld")
faces = data["faces"]

con = ConnectivityDS(faces,10)

verticies = []
for i in VertexVRing(1,con)
    push!(verticies,i)
end

@test sort(verticies)==[2,4,5,14,15]

triverticies = []
for i in DoubleVertexVRing(1,con)
    push!(triverticies,i)
end

@test (2,14) in triverticies

# Here I will also have tests for connectivity table

##### Surface Properties ######

@info "Tests for surface properties"

using JLD
data = load(joinpath(dirname(@__FILE__),"sphere.jld"))
#data = load("sphere.jld")
points = data["points"]
faces = data["faces"]

n = Array{Float64}(undef,size(points)...)
NormalVectors!(n,points,faces,i->FaceVRing(i,faces))
curvatures = Array{Float64}(undef,size(points,2))
MeanCurvatures!(curvatures,points,faces,n,i->FaceVRing(i,faces))

angle_ = 0
curvaturer = 0

for xkey in 1:size(points,2)
    iter = FaceVRing(xkey,faces)
    ncalc = n[:,xkey]
    ccalc = curvatures[xkey]

    nT = points[:,xkey]/norm(points[:,xkey])
    angl = acos(dot(ncalc,nT))
    if angl>angle_
        global angle_ = angl
    end

    cT = 1.
    N = size(points,2)
    global curvaturer += abs((ccalc - cT)/cT)/N
end

@test angle_*180/pi < 1.3
@test curvaturer < 0.01

@test isapprox(volume(points,faces),4/3*pi,atol=0.05)
@test isapprox(*(FitEllipsoid(points)...),1)

### And now integrators
@info "Testing pasive mesh stabilisation"

function velocity(t,pos)
    x,y,z = pos

    x = x*0.15 + 0.35
    y = y*0.15 + 0.35
    z = z*0.15 + 0.35

    u = 2*sin(pi*x)^2 * sin(2*pi*y) * sin(2*pi*z) * sin(2/3*pi*t)
    v = - sin(2*pi*x) * sin(pi*y)^2 * sin(2*pi*z) * sin(2/3*pi*t)
    w = - sin(2*pi*x) * sin(2*pi*y) * sin(pi*z)^2 * sin(2/3*pi*t)

    [u,v,w] /0.15
end

v = zero(points)
t = 1
for i in 1:size(points,2)
    v[:,i] = velocity(t,points[:,i])
end

res =  stabilise(points,faces,v)
println("Energy before minimization Finit=$(res.Finit) after Fres=$(res.Fres)")


text:
"""
    CorrelatedNoiseSource(eltype, samplerate, nchannels, std, correlation)

CorrelatedNoiseSource is a two-channel  noise signal generator with controlled correlation between channels.

Inputs
------
* `samplerate` specifies the sample rate of the signal.
* `nchannels` specifies the number of channels of the signal.
* `std` specifies the desired standard deviation of the signal.
* `correlation` specifies the desired correlation between the signals.

Output
------
* SampleSource object

Example
-------
```julia
source_object = CorrelatedNoiseSource(Float64, 48000, 2, 0.3, 0.75)
cn = read(source_object, 480)         # Specify number of samples of signal to generate
cn = read(source_object, 50u"ms")     # Specify length of time of signal to generate
```
"""
mutable struct CorrelatedNoiseSource{T} <: SampleSource
    samplerate::Float64
    nchannels::Int64
    cholcov::Array{Float64}

    function CorrelatedNoiseSource(eltype, samplerate::Number, nchannels::Number, std::Number, corr::Number)

        @assert nchannels==2 "Only two channels are supported for CorrelatedNoiseSource"

        if corr < 1
            correlation_matrix  = [1.0 corr ;  corr 1.0]
            standard_deviation  = [std 0.0 ;  0.0 std]
            covariance_matrix   = standard_deviation*correlation_matrix*standard_deviation
            cholcov = cholesky(covariance_matrix)
            cholcov = cholcov.U
        else
            cholcov = [std std ; 0 0]
        end

        new{eltype}(Float64(samplerate), Int64(nchannels), Array{Float64}(cholcov))
    end

    CorrelatedNoiseSource(eltype, samplerate::Unitful.Frequency, nchannels::Number, std::Number, corr::Number) = CorrelatedNoiseSource(eltype, samplerate |> u"Hz" |> ustrip, nchannels, std, corr)
end

Base.eltype(::CorrelatedNoiseSource{T}) where T = T
nchannels(source::CorrelatedNoiseSource) = source.nchannels
samplerate(source::CorrelatedNoiseSource) = source.samplerate

function unsafe_read!(source::CorrelatedNoiseSource, buf::Array, frameoffset, framecount)
    buf[1+frameoffset:framecount+frameoffset, 1:source.nchannels] = randn(framecount, source.nchannels) * source.cholcov
    framecount
end


text:
using VariableSelector
using Documenter

DocMeta.setdocmeta!(VariableSelector, :DocTestSetup, :(using VariableSelector); recursive=true)

makedocs(;
    modules=[VariableSelector],
    authors="rbac <ronbactawar@gmail.com> and contributors",
    repo="https://github.com/rbac/VariableSelector.jl/blob/{commit}{path}#{line}",
    sitename="VariableSelector.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rbac.github.io/VariableSelector.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rbac/VariableSelector.jl",
)

max_stars_repo_name: ronbactawar/VariableSelection

max_issues_repo_name: ronbactawar/VariableSelection

max_forks_repo_name: ronbactawar/VariableSelection


text:
module OpenML

using HTTP
using JSON
import ARFFFiles
using Markdown
if VERSION > v"1.3.0"
    using Pkg.Artifacts
end

export OpenML

include("data.jl")

end # module


text:
module RocketAccumulatedOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: accumulated()" begin

    println("Testing: operator accumulated()")

    run_proxyshowcheck("Accumulated", accumulated())

    run_testset([
        (
            source      = from(1:3) |> accumulated(),
            values      = @ts([ [ 1 ], [ 1, 2 ], [ 1, 2, 3 ], c ]),
            source_type = Vector{Int}
        ),
        (
            source      = from([ "a", "b", "c" ]) |> accumulated(),
            values      = @ts([ [ "a" ], [ "a", "b" ], [ "a", "b", "c" ], c ]),
            source_type = Vector{String}
        ),
        (
            source      = completed(Int) |> accumulated(),
            values      = @ts(c),
            source_type = Vector{Int}
        ),
        (
            source      = faulted(Int, "e") |> accumulated(),
            values      = @ts(e("e")),
            source_type = Vector{Int}
        ),
        (
            source      = never(Int) |> accumulated(),
            values      = @ts(),
            source_type = Vector{Int}
        )
    ])

end

end


text:
abstract type MultiFieldStyle end

struct ConsecutiveMultiFieldStyle <: MultiFieldStyle end

struct StridedMultiFieldStyle <: MultiFieldStyle end

"""
    struct MultiFieldFESpace{S<:MultiFieldStyle,B} <: FESpace
      spaces::Vector{<:SingleFieldFESpace}
      multi_field_style::S
      constraint_style::Val{B}
    end
"""
struct MultiFieldFESpace{MS<:MultiFieldStyle,CS<:ConstraintStyle,V} <: FESpace
  vector_type::Type{V}
  spaces::Vector{<:SingleFieldFESpace}
  multi_field_style::MS
  constraint_style::CS
  function MultiFieldFESpace(
    ::Type{V},
    spaces::Vector{<:SingleFieldFESpace},
    multi_field_style::MultiFieldStyle) where V
    @assert length(spaces) > 0

    MS = typeof(multi_field_style)
    if any( map(has_constraints,spaces) )
      constraint_style = Constrained()
    else
      constraint_style = UnConstrained()
    end
    CS = typeof(constraint_style)
    new{MS,CS,V}(V,spaces,multi_field_style,constraint_style)
  end
end

"""
    MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
"""
function MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
  Ts = map(get_dof_value_type,spaces)
  T = typeof(*(map(zero,Ts)...))
  MultiFieldFESpace(Vector{T},spaces,ConsecutiveMultiFieldStyle())
end

function MultiFieldFESpace(::Type{V},spaces::Vector{<:SingleFieldFESpace}) where V
  MultiFieldFESpace(V,spaces,ConsecutiveMultiFieldStyle())
end

MultiFieldStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = S()
MultiFieldStyle(f::MultiFieldFESpace) = MultiFieldStyle(typeof(f))

# Implementation of FESpace

function FESpaces.get_triangulation(f::MultiFieldFESpace)
  s1 = first(f.spaces)
  trian = get_triangulation(s1)
  @check all(map(i->trian===get_triangulation(i),f.spaces))
  trian
end

function FESpaces.num_free_dofs(f::MultiFieldFESpace)
  n = 0
  for U in f.spaces
    n += num_free_dofs(U)
  end
  n
end

FESpaces.get_dof_value_type(f::MultiFieldFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V)

FESpaces.get_vector_type(f::MultiFieldFESpace) = f.vector_type

FESpaces.ConstraintStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = B()

function FESpaces.get_cell_shapefuns(f::MultiFieldFESpace)

  # Compute the cell axes
  blocks = map(get_cell_shapefuns,f.spaces)
  blocks_data = map(get_data,blocks)
  cell_axes_blocks = map(i->lazy_map(axes,i),blocks_data)
  cell_axes = lazy_map(_multifield_axes_dofs,cell_axes_blocks...)

  # Compute the underlying single-fields
  nfields = length(f.spaces)
  all_bases = map(1:nfields) do i
    bsize = (nfields,)
    dv = blocks[i]
    cell_basis = lazy_map(Fields.BlockFieldArrayCooMap(bsize,[(i,)]),cell_axes,get_data(dv))
    FEBasis(cell_basis,get_triangulation(dv),TestBasis(),DomainStyle(dv))
  end
  MultiFieldCellField(all_bases)
end

function _multifield_axes_dofs(axs...)
  rs = map(first,axs)
  (append_ranges([rs...]),)
end

function FESpaces.get_cell_shapefuns_trial(f::MultiFieldFESpace)
  blocks = get_cell_shapefuns(f).single_fields
  nfields = length(blocks)
  all_bases = map(1:nfields) do i
    dv = blocks[i]
    cell_basis = lazy_map(transpose,get_data(dv))
    FEBasis(cell_basis,get_triangulation(dv),TrialBasis(),DomainStyle(dv))
  end
  MultiFieldCellField(all_bases)
end

function FESpaces.FEFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    FEFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function FESpaces.EvaluationFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    EvaluationFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function CellData.CellField(fe::MultiFieldFESpace,cell_values)
  single_fields = map(1:length(fe.spaces)) do i
    cell_values_field = lazy_map(a->a[Block(i)],cell_values)
    CellField(fe.spaces[i],cell_values_field)
  end
  MultiFieldCellField(single_fields)
end

function CellData.CellField(
  fe::MultiFieldFESpace,cell_values::LazyArray{<:Fill{<:BlockArrayCooMap}})
  single_fields = map(1:length(fe.spaces)) do i
    ma = cell_values.maps.value
    cell_values_field = Fields._get_cell_block(ma,cell_values,Block(i))
    CellField(fe.spaces[i],cell_values_field)
  end
  MultiFieldCellField(single_fields)
end

"""
    restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
"""
function restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
  _restrict_to_field(f,MultiFieldStyle(f),free_values,field)
end

function  _restrict_to_field(f,::MultiFieldStyle,free_values,field)
  @notimplemented
end

function  _restrict_to_field(f,::ConsecutiveMultiFieldStyle,free_values,field)
  offsets = compute_field_offsets(f)
  U = f.spaces
  pini = offsets[field] + 1
  pend = offsets[field] + num_free_dofs(U[field])
  SubVector(free_values,pini,pend)
end

"""
    compute_field_offsets(f::MultiFieldFESpace)
"""
function compute_field_offsets(f::MultiFieldFESpace)
  @assert MultiFieldStyle(f) == ConsecutiveMultiFieldStyle()
  U = f.spaces
  n = length(U)
  offsets = zeros(Int,n)
  for i in 1:(n-1)
    Ui = U[i]
    offsets[i+1] = offsets[i] + num_free_dofs(Ui)
  end
  offsets
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace)
  data = map(get_cell_isconstrained,f.spaces)
  lazy_map( (args...) -> +(args...)>0,  data...)
end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace)
  blocks = map(get_cell_constraints,f.spaces)
  blockids = [ (i,i) for i in 1:length(f.spaces)]
  nfields = length(f.spaces)
  bsize = (nfields,nfields)

  blocks_axes = map( i->lazy_map(axes,i), blocks)

  function multifield_axes_constraints(axs...)
    rs1 = map(j->j[1],axs)
    rs2 = map(j->j[2],axs)
    r1 = append_ranges([rs1...])
    r2 = append_ranges([rs2...])
    (r1,r2)
  end
  cell_axes = lazy_map(multifield_axes_constraints,blocks_axes...)
  lazy_map(BlockArrayCooMap(bsize,blockids),cell_axes,blocks...)
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace)
  get_cell_dof_ids(f,MultiFieldStyle(f))
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,::MultiFieldStyle)
  @notimplemented
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,::ConsecutiveMultiFieldStyle)
  offsets = compute_field_offsets(f)
  spaces = f.spaces
  function fun(i,space)
    cell_dofs = get_cell_dof_ids(space)
    if i == 1
      return cell_dofs
    end
    offset = Int32(offsets[i])
    o = Fill(offset,length(cell_dofs))
    lazy_map(Broadcasting(_sum_if_first_positive),cell_dofs,o)
  end
  blocks = [ fun(i,space) for (i,space) in enumerate(spaces) ]
  bids = [ (i,) for i in 1:length(spaces)]
  bsize = (length(f.spaces),)
  cell_axes_blocks = map(i->lazy_map(axes,i),blocks)
  cell_axes = lazy_map(_multifield_axes_dofs,cell_axes_blocks...)
  lazy_map(BlockArrayCooMap(bsize,bids),cell_axes,blocks...)
end

function _sum_if_first_positive(a,b)
  if a>0
    return a+b
  else
    return a
  end
end

# API for multi field case

"""
    num_fields(f::MultiFieldFESpace)
"""
function num_fields(f::MultiFieldFESpace)
  length(f.spaces)
end

Base.iterate(m::MultiFieldFESpace) = iterate(m.spaces)

Base.iterate(m::MultiFieldFESpace,state) = iterate(m.spaces,state)

Base.getindex(m::MultiFieldFESpace,field_id::Integer) = m.spaces[field_id]

# API for the ConsecutiveMultiFieldStyle
import Gridap.FESpaces: interpolate
import Gridap.FESpaces: interpolate_everywhere
import Gridap.FESpaces: interpolate_dirichlet

@deprecate(
  interpolate(fs::MultiFieldFESpace, object),
  interpolate(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_everywhere(fs::MultiFieldFESpace, object),
  interpolate_everywhere(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_dirichlet(fs::MultiFieldFESpace, object),
  interpolate_dirichlet(object, fs::MultiFieldFESpace)
)

"""
The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs
even in the case that the given cell field does not fulfill them)
"""
function FESpaces.interpolate(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    uhi = interpolate!(object, free_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
like interpolate, but also compute new degrees of freedom for the dirichlet component.
The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
"""
function FESpaces.interpolate_everywhere(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_everywhere!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
"""
function FESpaces.interpolate_dirichlet(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_dirichlet!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end


text:
using Flux: testmode!

@testset "Dropout" begin
  x = [1.,2.,3.]
  @test x == testmode!(Dropout(0.1))(x)
  @test x == Dropout(0)(x)
  @test zeros(x) == Dropout(1)(x)

  x = rand(100)
  m = Dropout(0.9)
  y = m(x)
  @test count(a->a==0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a==0, y) == 0
  testmode!(m, false)
  y = m(x)
  @test count(a->a==0, y) > 50

  x = rand(100)
  m = Chain(Dense(100,100),
            Dropout(0.9))
  y = m(x)
  @test count(a->a == 0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a == 0, y) == 0
end


text:
using FFTW
using BenchmarkTools
using HDF5
using Tullio

itot = 256; jtot = 192; ktot = 128;

float_type = Float32

x = collect(float_type, 0:itot-1) ./ itot .* 2pi
y = collect(float_type, 0:jtot-1) ./ jtot .* 2pi
z = collect(float_type, 0:ktot-1) ./ ktot .* 2pi

a = zeros(float_type, itot, jtot, ktot)
a_new = zeros(float_type, itot, jtot, ktot)
tmp = zeros(float_type, itot, jtot, ktot)

@tullio a[i, j, k] = cos(4x[i]) + cos(8y[j])

FFTW.set_num_threads(Threads.nthreads())

fft_plan_f = FFTW.plan_r2r(a, FFTW.R2HC, (1, 2), flags=FFTW.ESTIMATE)
fft_plan_b = FFTW.plan_r2r(a, FFTW.HC2R, (1, 2), flags=FFTW.ESTIMATE)

function forward_backward!(tmp, a_new, a, fft_plan_f, fft_plan_b, itot, jtot)
    tmp[:, :, :] = fft_plan_f * a
    a_new[:, :, :] = fft_plan_b * tmp ./ (itot * jtot)
end

@btime forward_backward!(
    $tmp, $a_new, $a, $fft_plan_f, $fft_plan_b, $itot, $jtot)
forward_backward!(tmp, a_new, a, fft_plan_f, fft_plan_b, itot, jtot)

h5open("fft_test.h5", "w") do fid
    write(fid, "tmp", tmp[:, :, :])
    write(fid, "a", a[:, :, :])
    write(fid, "a_new", a_new[:, :, :])
end


text:
using DataDeps

register(DataDep(
    "vowel",
    """
    Dataset: vowel
    Website: https://web.stanford.edu/~hastie/ElemStatLearn/

    Observations: 990 (528 + 462)
    Features:     10
    Classes:      11

    This info is the original source information for these data.

        NAME: Vowel Recognition (Deterding data)

        SUMMARY: Speaker independent recognition of the eleven steady state vowels
        of British English using a specified training set of lpc derived log area
        ratios.

        SOURCE: David Deterding  (data and non-connectionist analysis)
                Mahesan Niranjan (first connectionist analysis)
                Tony Robinson    (description, program, data, and results)

        To contact Tony Robinson by electronic mail, use address
        "ajr@dsl.eng.cam.ac.uk"

        MAINTAINER: neural-bench@cs.cmu.edu

    See the complete information at https://web.stanford.edu/~hastie/ElemStatLearn/datasets/vowel.info.txt
    """,
    [
    "https://web.stanford.edu/~hastie/ElemStatLearn/datasets/vowel.train",
    "https://web.stanford.edu/~hastie/ElemStatLearn/datasets/vowel.test"
    ],
    "f6f7d76de6247e85add6ec3e0d5988f7406338f1562897e558844ae750b6b6e0",
    #
    # Disclaimer: Want to combine both .train and .test into a single .csv.
    # Design of post_fetch_method may not have been designed with this in mind.
    # Workaround here is to do nothing to first file, then assume first file has been
    # downloaded and can be loaded in the second function.
    # This works in practice, but will fail when running preupload_check.
    #
    # Reason: `path` is randomly generated and differs between the two calls.
    #
    post_fetch_method = [
    identity,
    path -> begin
        # Read both training and testing files.
        dir = dirname(path)
        tra = joinpath(dir, "vowel.train")
        tes = joinpath(dir, "vowel.test")
        df = vcat(
            CSV.read(tra, DataFrame, header=true),
            CSV.read(tes, DataFrame, header=true)
        )

        # Process the data as usual.
        MVDA.process_dataset(df,
            target_index=2,
            feature_indices=3:ncol(df),
            ext=".csv",
        )

        # Clean up by removing separate training and testing files.
        rm(tra)
        rm(tes)
    end],
))


text:
using MPThemes
using Test

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


text:
using Test
using Codex.TransformExport
using Codex.Questionnaires

@testset "Questionnaires" begin
    include("personality.jl")

    data_dir = joinpath(Codex.PROJECT_ROOT, "test", "data", "2020-08")

    df = Codex.Questionnaires.responses(data_dir, "second")
    @test size(df) == (3, 8)

    df = Codex.Questionnaires.responses(data_dir, "delta")
    @test df[1, :score] == 26.0

    df = Codex.Questionnaires.responses(data_dir, "foxtrot")
    @test df[:, [:correct, :incorrect, :dontknow]] ==
        DataFrame(correct = [8], incorrect = [3], dontknow = [0])

    df = Codex.Questionnaires.responses(data_dir, "kilo")
    @test df[:, [:optimism, :pessimism]] == DataFrame(optimism = [10], pessimism = [7])

    df = Codex.Questionnaires.responses(data_dir, "mike")
    @test df[1, :optimism] ≈ 19 / 6

    responses_dir = joinpath(data_dir, "responses")
    id = "aaab"
    unfinished = Codex.Questionnaires.unfinished_questionnaires(responses_dir, id)
    @test unfinished == ["first", "personality", "third"]

    ids = Codex.Questionnaires.all_ids(responses_dir)
    @test ids == ["aaaa", "aaab", "aaac", "aaad"]

    U = Codex.Questionnaires.unfinished_info(responses_dir; required=["first"])
    @test getproperty.(U, :id) == ["aaab", "aaac", "aaad"]
    @test getproperty.(U, :unfinished_questionnaires) == fill(["first"], 3)
end


text:
using Documenter, DynamicGrids, KernelAbstractions

makedocs(
    modules = [DynamicGrids],
    sitename = "DynamicGrids.jl",
    checkdocs = :all,
    strict = true,
)

deploydocs(
    repo = "github.com/cesaraustralia/DynamicGrids.jl.git",
)


text:
using BSON

function 사전가져오기(path::String=normpath(@__DIR__, "..", "data", "dict.bson"))
    BSON.@load path dict
    return dict
end


text:
module ReinforcementLearningEnvironmentClassicControl
using GR
using ReinforcementLearningBase
import ReinforcementLearningBase: interact!, getstate, reset!, actionspace, plotenv

include("cartpole.jl")
include("mountaincar.jl")
include("pendulum.jl")

export CartPole, MountainCar, Pendulum

end # module


text:
module LightXML

    export

    # common
    name, free,

    # nodes
    AbstractXMLNode,
    XMLAttr, XMLAttrIter, XMLNode, XMLNodeIter, XMLElement, XMLElementIter,
    nodetype, value, content, attribute, has_attribute,
    is_elementnode, is_textnode, is_commentnode, is_cdatanode, is_blanknode,
    child_nodes, has_children, attributes, has_attributes, attributes_dict,
    child_elements, find_element, get_elements_by_tagname,

    new_element, add_child, new_child, new_textnode, add_text, add_cdata,
    set_attribute, set_attributes,

    # document
    XMLDocument, version, encoding, compression, standalone, root,
    parse_file, parse_string, save_file, set_root, create_root

    include("clib.jl")
    include("errors.jl")

    include("utils.jl")
    include("nodes.jl")
    include("document.jl")
    include("cdata.jl")
end


text:
################################################################################
#
#   Calculates integer coordinates in ℤⁿ lying inside the n
#          dimensional ellipsoid
#
#   Implementation of Remarks page 1426 of
#   B. Deconinck, M.  Heil, A. Bobenko, M. van Hoeij and M. Schmies,
#    Computing Riemann Theta Functions, Mathematics of Computation, 73, (2004),
#    1417-1442.
#
################################################################################

const padding = 0.5 # necessary for uniform approximation (page 1428)

function innerpoints(T::Matrix{Float64}, radius::Float64)
    n   = size(T, 1) # problem size

    # precalculate what is independent of points
    δcs = [ inv(T[1:i-1,1:i-1]) * T[1:i-1,i] for i in 1:n ]
    Tgg = diag(T)

    ns     = Vector{Float64}(undef, n)
    points = Vector{Float64}[]
    function _innerpoints(Rₒ, c, g)
        hw = Rₒ / Tgg[g]
        for ng in ceil(c[g]-hw):floor(c[g]+hw)
            ns[g] = ng
            if g == 1
                push!(points, copy(ns))
            else
                δcn = (ng - c[g])
                nc = c[1:g-1] - δcs[g] * δcn
                nsqRₒ = Rₒ^2 - (Tgg[g] * δcn)^2
                nsqRₒ > 0. && _innerpoints(sqrt(nsqRₒ), nc, g-1)
            end
        end
    end

    # add to initial R such that half-width 'hw' is augmented by 'padding'
    padded_Rₒ = radius / sqrt(π) + padding * Tgg[n]
    _innerpoints(padded_Rₒ, zeros(n), n)

    points
end


text: # should we try to install lme4 or even afex automatically at build time?


text:
module TikzCDs

using TikzPictures

import TikzPictures: PDF, TEX, TIKZ, SVG, save, TikzDocument, push!, tikzCommand
export PDF, TEX, TIKZ, SVG, save, TikzDocument, push!, tikzCommand
import LaTeXStrings: LaTeXString, @L_str, @L_mstr
export LaTeXString, @L_str, @L_mstr

export TikzCD

function TikzCD(data::AbstractString; options="", preamble="", enableWrite18=true)
  # check for tikz-cd package being imported
  if !occursin("\\usepackage{tikz-cd}", preamble)
    preamble = strip(string("\\usepackage{tikz-cd}\n", preamble))
  end
  # tikz-cd doesn't support the content being in a $$, so this helps when using LaTeXStrings
  if typeof(data) == LaTeXString
    data = strip(data[2:end-1])
  end
  TikzPicture(data, options=options, preamble=preamble, environment="tikzcd", enableWrite18=enableWrite18)
end

end


text:
using Test
using PlayingCards52

@test true

for k = 1:52
    C = Card(k)
    kk = index(C)
    @test k == kk
end

ace_high(false)
x = deck()
sort!(x)
@test index(x[1]) == 1

AC = Card(:clubs, 1)
twoD = Card(:diamonds, 2)

ace_high(true)
@test twoD < AC
ace_high(false)
@test twoD > AC

@test color(Card(:clubs, 5)) == color(Card(:spades, 2))
@test color(Card(:clubs, 5)) != color(Card(:diamonds, 2))
@test color(Suit(1)) == color(♠)

@test Char(Card(:clubs, 5)) == '🃕'

d = collect(1:10)
cut!(d, 4)
cut!(d, 6)
@test d == collect(1:10)

riffle!(d)
@test sort(d) == collect(1:10)

# test alternative notation

@test A♣ == Card(:clubs, 1)
@test K♣ == 13♣
@test 7♠ == Card(:spades, 7)

@test 5♢ == 5♦
@test A♡ == A♥

@test Suit(:clubs) == ♣
@test ♡ == ♥

@test name(Card(:clubs, 13)) == "king of clubs"


text:
## --- ArrayStats.jl

## --- Filtering
    A = 1:100
    @test A[inpctile(A,80)] == 11:90

## --- NaN handling functions, simple cases
    A = [1:10; fill(NaN,10)]
    B = [fill(NaN,10); 11:20]
    @test countnans(A) == countnotnans(A) == 10
    @test A[nanmask(A)] == 1:10
    @test nanadd(A,B) == 1:20
    @test nanadd!(A,B) == 1:20
    @test zeronan!(B) == [fill(0,10); 11:20]

## --- Summary statistics: simple cases, Float64
    A = [1:10.0..., NaN]
    @test nansum(A) == 55.0
    @test nanmean(A) == 5.5
    @test nanmean(A, ones(11)) == 5.5 # weighted
    @test nanrange(A) == 9.0
    @test nanminimum(A) == 1.0
    @test nanmaximum(A) == 10.0
    @test nanextrema(A) == (1.0, 10.0)
    @test nanvar([1,2,3,NaN]) == 1.0
    @test nanstd([1,2,3,NaN]) == 1.0
    @test nanstd([1,2,3,NaN], ones(4)) == 1.0 # weighted
    @test nanmad([1,2,3,NaN]) == 1.0
    @test nanaad([1,2,3,NaN]) ≈ 2/3
    @test nanmedian([1,2,3,NaN]) == 2.0
    @test nanpctile([0:100...,NaN],99) == 99.0
    @test nanmin(1.,2.) == 1.
    @test nanmax(1.,2.) == 2.

## --- Arrays containing only NaNs should yield NaN
    A = fill(NaN,10)
    @test nansum(A) == 0
    @test isnan(nanmean(A))
    @test isnan(nanmean(A, ones(10))) # weighted
    @test isnan(nanrange(A))
    @test isnan(nanminimum(A))
    @test isnan(nanmaximum(A))
    @test all(isnan.(nanextrema(A)))
    @test isnan(nanvar(A))
    @test isnan(nanvar(A, mean=NaN))
    @test isnan(nanstd(A))
    @test isnan(nanstd(A, ones(10))) # weighted
    @test isnan(nanaad(A))
    @test isnan(nanmin(NaN,NaN))
    @test isnan(nanmax(NaN,NaN))

## --- Summary statistics: simple cases, Int64
    A = collect(1:10)
    @test nansum(A) == 55.0
    @test nanmean(A) == 5.5
    @test nanmean(A, ones(10)) == 5.5 # weighted
    @test nanrange(A) == 9.0
    @test nanminimum(A) == 1.0
    @test nanmaximum(A) == 10.0
    @test nanextrema(A) == (1.0, 10.0)
    @test nanvar([1,2,3]) == 1.0
    @test nanvar([1,2,3], mean=2) == 1.0
    @test nanstd([1,2,3]) == 1.0
    @test nanstd([1,2,3], ones(3)) == 1.0 # weighted
    @test nanmad([1,2,3]) == 1.0
    @test nanaad([1,2,3]) ≈ 2/3
    @test nanmedian([1,2,3]) == 2.0
    @test nanpctile([0:100...],99) == 99.0
    @test nanmin(1,2) == 1
    @test nanmax(1,2) == 2

## --- Summary statistics: simple cases, ranges
    A = 1:10
    @test nansum(A) == 55.0
    @test nanmean(A) == 5.5
    @test nanmean(A, ones(10)) == 5.5 # weighted
    @test nanrange(A) == 9.0
    @test nanminimum(A) == 1.0
    @test nanmaximum(A) == 10.0
    @test nanextrema(A) == (1.0, 10.0)
    @test nanvar(1:3) == 1.0
    @test nanvar(1:3, mean=2) == 1.0
    @test nanstd(1:3) == 1.0
    @test nanstd(1:3, ones(3)) == 1.0 # weighted
    @test nanmad(1:3) == 1.0
    @test nanaad(1:3) ≈ 2/3
    @test nanmedian(1:3) == 2.0
    @test nanpctile(0:100,99) == 99.0

    A = 1:10.
    @test nansum(A) == 55.0
    @test nanmean(A) == 5.5
    @test nanmean(A, ones(10)) == 5.5 # weighted
    @test nanrange(A) == 9.0
    @test nanminimum(A) == 1.0
    @test nanmaximum(A) == 10.0
    @test nanextrema(A) == (1.0, 10.0)
    @test nanvar(1:3.) == 1.0
    @test nanvar(1:3., mean=2.0) == 1.0
    @test nanstd(1:3.) == 1.0
    @test nanstd(1:3., ones(3)) == 1.0 # weighted
    @test nanmad(1:3.) == 1.0
    @test nanaad(1:3.) ≈ 2/3
    @test nanmedian(1:3.) == 2.0
    @test nanpctile(0:100.,99) == 99.0

## --- Summary statistics: dimensional tests, Int64
    A = collect(reshape(1:300,100,3))
    @test nansum(A, dims=1) == sum(A, dims=1)
    @test nansum(A, dims=2) == sum(A, dims=2)
    @test nanminimum(A, dims=1) == minimum(A, dims=1)
    @test nanminimum(A, dims=2) == minimum(A, dims=2)
    @test nanmaximum(A, dims=1) == maximum(A, dims=1)
    @test nanmaximum(A, dims=2) == maximum(A, dims=2)
    @test nanextrema(A, dims=1) == [(1.0, 100.0) (101.0, 200.0) (201.0, 300.0)]
    @test nanmean(A, dims=1) == mean(A, dims=1)
    @test nanmean(A, dims=2) == mean(A, dims=2)
    @test nanmean(A, ones(size(A)), dims=1) == mean(A, dims=1) # weighted
    @test nanmean(A, ones(size(A)), dims=2) == mean(A, dims=2) # weighted
    @test nanvar(A, dims=1) ≈ var(A, dims=1)
    @test nanvar(A, dims=2) ≈ var(A, dims=2)
    @test nanstd(A, dims=1) ≈ std(A, dims=1)
    @test nanstd(A, dims=2) ≈ std(A, dims=2)
    @test nanstd(A, dims=1, mean=nanmean(A,dims=1)) ≈ std(A, dims=1)
    @test nanstd(A, dims=2, mean=nanmean(A,dims=2)) ≈ std(A, dims=2)
    @test nanstd(A, ones(size(A)), dims=1) ≈ std(A, dims=1) # weighted
    @test nanstd(A, ones(size(A)), dims=2) ≈ std(A, dims=2) # weighted
    @test nanmad(A, dims=1) == [25.0 25.0 25.0]
    @test nanmad(A, dims=2) == fill(100.0,100,1)
    @test nanaad(A, dims=1) == [25.0 25.0 25.0] #
    @test nanaad(A, dims=2) ≈ fill(200/3,100,1)
    @test nanmedian(A, dims=1) == median(A, dims=1)
    @test nanmedian(A, dims=2) == median(A, dims=2)
    @test nanpctile(A, 10, dims=1) ≈ [10.9  110.9  210.9]
    @test nanpctile(A, 10, dims=2) ≈ 21:120

## --- Summary statistics: dimensional tests, Float64
    A = collect(reshape(1:300.,100,3))
    @test nansum(A, dims=1) == sum(A, dims=1)
    @test nansum(A, dims=2) == sum(A, dims=2)
    @test nanminimum(A, dims=1) == minimum(A, dims=1)
    @test nanminimum(A, dims=2) == minimum(A, dims=2)
    @test nanmaximum(A, dims=1) == maximum(A, dims=1)
    @test nanmaximum(A, dims=2) == maximum(A, dims=2)
    @test nanextrema(A, dims=1) == [(1.0, 100.0) (101.0, 200.0) (201.0, 300.0)]
    @test nanmean(A, dims=1) == mean(A, dims=1)
    @test nanmean(A, dims=2) == mean(A, dims=2)
    @test nanmean(A, ones(size(A)), dims=1) == mean(A, dims=1) # weighted
    @test nanmean(A, ones(size(A)), dims=2) == mean(A, dims=2) # weighted
    @test nanvar(A, dims=1) ≈ var(A, dims=1)
    @test nanvar(A, dims=2) ≈ var(A, dims=2)
    @test nanstd(A, dims=1) ≈ std(A, dims=1)
    @test nanstd(A, dims=2) ≈ std(A, dims=2)
    @test nanstd(A, dims=1, mean=nanmean(A,dims=1)) ≈ std(A, dims=1)
    @test nanstd(A, dims=2, mean=nanmean(A,dims=2)) ≈ std(A, dims=2)
    @test nanstd(A, ones(size(A)), dims=1) ≈ std(A, dims=1) # weighted
    @test nanstd(A, ones(size(A)), dims=2) ≈ std(A, dims=2) # weighted
    @test nanmedian(A, dims=1) == median(A, dims=1)
    @test nanmedian(A, dims=2) == median(A, dims=2)
    @test nanminimum(A, dims=1) == [1 101 201]
    @test dropdims(nanminimum(A, dims=2), dims=2) == 1:100
    @test nanmaximum(A, dims=1) == [100 200 300]
    @test dropdims(nanmaximum(A, dims=2), dims=2) == 201:300
    @test nanmean(A, dims=1) == [50.5 150.5 250.5]
    @test dropdims(nanmean(A, dims=2), dims=2) == 101:200
    @test nanmean(A, ones(size(A)), dims=1) == [50.5 150.5 250.5] # weighted
    @test dropdims(nanmean(A, ones(size(A)), dims=2), dims=2) == 101:200 # weighted
    @test nanstd(A, dims=1) ≈ fill(29.011491975882016, 1, 3)
    @test nanstd(A, dims=2) ≈ fill(100, 100, 1)
    @test nanstd(A, ones(size(A)), dims=1) ≈ fill(29.011491975882016, 1, 3) # weighted
    @test nanstd(A, ones(size(A)), dims=2) ≈ fill(100, 100, 1) # weighted
    @test nanmedian(A, dims=1) == [50.5 150.5 250.5]
    @test dropdims(nanmedian(A, dims=2), dims=2) == 101:200
    @test nanpctile(A, 10, dims=1) ≈ [10.9 110.9 210.9]
    @test nanpctile(A, 10, dims=2) ≈ 21:120
    @test nanmad(A, dims=1) == [25.0 25.0 25.0]
    @test nanmad(A, dims=2) == fill(100.0, 100, 1)
    @test nanaad(A, dims=1) == [25.0 25.0 25.0]
    @test nanaad(A, dims=2) ≈ fill(200/3, 100, 1)

## --- Summary statistics: dimensional tests, Float64
    A = collect(reshape(1:300.,100,3))
    @test nansum(A, dim=1) == vec(sum(A, dims=1))
    @test nansum(A, dim=2) == vec(sum(A, dims=2))
    @test nanminimum(A, dim=1) == vec(minimum(A, dims=1))
    @test nanminimum(A, dim=2) == vec(minimum(A, dims=2))
    @test nanmaximum(A, dim=1) == vec(maximum(A, dims=1))
    @test nanmaximum(A, dim=2) == vec(maximum(A, dims=2))
    @test nanextrema(A, dim=1) == [(1.0, 100.0), (101.0, 200.0), (201.0, 300.0)]
    @test nanmean(A, dim=1) == vec(mean(A, dims=1))
    @test nanmean(A, dim=2) == vec(mean(A, dims=2))
    @test nanmean(A, ones(size(A)), dim=1) == vec(mean(A, dims=1)) # weighted
    @test nanmean(A, ones(size(A)), dim=2) == vec(mean(A, dims=2)) # weighted
    @test nanvar(A, dim=1) ≈ vec(var(A, dims=1))
    @test nanvar(A, dim=2) ≈ vec(var(A, dims=2))
    @test nanvar(A, dim=1, mean=nanmean(A,dims=1)) ≈ vec(var(A, dims=1))
    @test nanvar(A, dim=2, mean=nanmean(A,dims=2)) ≈ vec(var(A, dims=2))
    @test nanstd(A, dim=1) ≈ vec(std(A, dims=1))
    @test nanstd(A, dim=2) ≈ vec(std(A, dims=2))
    @test nanstd(A, ones(size(A)), dim=1) ≈ vec(std(A, dims=1)) # weighted
    @test nanstd(A, ones(size(A)), dim=2) ≈ vec(std(A, dims=2)) # weighted
    @test nanmedian(A, dim=1) == vec(median(A, dims=1))
    @test nanmedian(A, dim=2) == vec(median(A, dims=2))
    @test nanminimum(A, dim=1) == [1, 101, 201]
    @test nanminimum(A, dim=2) == 1:100
    @test nanmaximum(A, dim=1) == [100, 200, 300]
    @test nanmaximum(A, dim=2) == 201:300
    @test nanmean(A, dim=1) == [50.5, 150.5, 250.5]
    @test nanmean(A, dim=2) == 101:200
    @test nanmean(A, ones(size(A)), dim=1) == [50.5, 150.5, 250.5] # weighted
    @test nanmean(A, ones(size(A)), dim=2) == 101:200 # weighted
    @test nanstd(A, dim=1) ≈ fill(29.011491975882016, 3)
    @test nanstd(A, dim=2) ≈ fill(100, 100)
    @test nanstd(A, ones(size(A)), dim=1) ≈ fill(29.011491975882016, 3) # weighted
    @test nanstd(A, ones(size(A)), dim=2) ≈ fill(100, 100) # weighted
    @test nanmedian(A, dim=1) == [50.5, 150.5, 250.5]
    @test nanmedian(A, dim=2) == 101:200
    @test nanpctile(A, 10, dim=1) ≈ [10.9, 110.9, 210.9]
    @test nanpctile(A, 10, dim=2) ≈ 21:120
    @test nanmad(A, dim=1) == [25.0, 25.0, 25.0]
    @test nanmad(A, dim=2) == fill(100.0, 100)
    @test nanaad(A, dim=1) == [25.0, 25.0, 25.0]
    @test nanaad(A, dim=2) ≈ fill(200/3, 100)

## --- Test fallbacks for complex reductions
    A = randn((2 .+ (1:6))...);
    @test nansum(A, dims=(4,5,6)) ≈ sum(A, dims=(4,5,6))
    @test nanmean(A, dims=(4,5,6)) ≈ mean(A, dims=(4,5,6))
    @test nanstd(A, dims=(4,5,6)) ≈ std(A, dims=(4,5,6))
    @test nanstd(A, dims=(4,5,6)) ≈ nanstd(A, dims=(4,5,6), mean=nanmean(A, dims=(4,5,6)))

## --- Test in-place vs. out-of-place versions

    A = rand(100)
    @test nanpctile(A, 25) == nanpctile!(A, 25)
    @test nanmedian(A) == nanmedian!(A)
    @test nanmad(A) == nanmad!(A)

## --- A few tests with other types
    for T in (Bool, Float16, Float32)
        let A = rand(T, 100)
            @test nansum(A) ≈ sum(A)
            @test nanmean(A) ≈ mean(A)
            @test nanstd(A) ≈ std(A)
        end
    end

    for T in (UInt8, UInt16, UInt32, Int8, Int16, Int32)
        let A = collect(T, 1:100)
            @test nansum(A) == 5050
            @test nanmean(A) == 50.5
            @test nanstd(A) ≈ 29.011491975882016
        end
    end

## --- Standardization

    @test nanstandardize!(collect(1:10.)) ≈ ((1:10) .- mean(1:10)) / std(1:10)
    @test nanstandardize(1:10.) ≈ ((1:10) .- mean(1:10)) / std(1:10)

## --- Moving averages

    # Moving average: 1D
    @test movmean(collect(1:10.),5) == movmean(1:10,5)
    @test movmean(1:10,5) == [2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.5, 9.0]
    # Moving average: 2D
    @test movmean(repeat(1:10,1,10),5) == repeat([2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.5, 9.0],1,10)
    @test movmean(repeat(1:10,1,10)',5) == repeat([2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.5, 9.0],1,10)'
    @test movmean([3 5 10 3; 4 2 5 8; 5 6 8 8; 2 6 10 6], 3) ≈ [7//2 29//6 11//2 13//2; 25//6 16//3 55//9 7//1; 25//6 16//3 59//9 15//2; 19//4 37//6 22//3 8//1]

## ---


text:
# Julia wrapper for header: /usr/local/include/libavutil/file.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_file_map,
    av_file_unmap,
    av_tempfile

function av_file_map(filename,bufptr,size,log_offset::Integer,log_ctx)
    ccall((:av_file_map,libavutil),Cint,(Cstring,Ptr{Ptr{UInt8}},Ptr{Csize_t},Cint,Ptr{Void}),filename,bufptr,size,log_offset,log_ctx)
end

function av_file_unmap(bufptr,size::Csize_t)
    ccall((:av_file_unmap,libavutil),Void,(Ptr{UInt8},Csize_t),bufptr,size)
end

function av_tempfile(prefix,filename,log_offset::Integer,log_ctx)
    ccall((:av_tempfile,libavutil),Cint,(Cstring,Ptr{Cstring},Cint,Ptr{Void}),prefix,filename,log_offset,log_ctx)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module dapr
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module proto
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    module sentry
      const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
      module v1
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("sentry_pb.jl")
      end
    end
  end
end


text:
module IntervalUnionArithmetic

using Reexport

@reexport using IntervalArithmetic

import Base: getindex, ∪, intersect, \, in, isequal
import Base: +, -, *, /, min, max, ^, log, <, >, exp, sin, cos, tan, sqrt
import IntervalArithmetic: hull, bisect, intersect, ⊆, ⊂

abstract type IntervalUnion{T} <: AbstractInterval{T} end

export
    IntervalUnion, IntervalU, intervalU, remove_empties, condense, closeGaps!,
    iscondensed, condense_weak, iscondensed_weak, env, hull, complement, left, right, bisect,
    intersect, setdiff, \, ⊆, ⊂, isequal, in

global MAXINTS = [100]

include("interval_unions.jl")
include("arithmetic.jl")
include("utilities.jl")
include("set_operations.jl")

end # module


text:
# Not used from v0.3.0 on (PowerModels v0.15.3) - Kept for sentimental reasons
function get_solution_acdc(pm::_PM.AbstractPowerModel, sol::Dict{String,Any})
    _PM.add_setpoint_bus_voltage!(sol, pm)
    _PM.add_setpoint_generator_power!(sol, pm)
    _PM.add_setpoint_branch_flow!(sol, pm)
    _PM.add_setpoint_dcline_flow!(sol, pm)
    if haskey(pm.data, "convdc")
        add_dc_bus_voltage_setpoint(sol, pm)
        add_dcconverter_setpoint(sol, pm)
        add_dcgrid_flow_setpoint(sol,pm)
        add_dcbranch_losses(sol, pm)
        add_dcconverter_losses(sol, pm)
        add_dcconverter_voltage_setpoint(sol, pm)
        add_dcconverter_firing_angle(sol, pm)
    end
    return sol
end

function add_dcconverter_setpoint(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc", "pgrid", :pconv_tf_fr)
    _PM.add_setpoint!(sol, pm, "convdc", "qgrid", :qconv_tf_fr)
    _PM.add_setpoint!(sol, pm, "convdc", "pconv", :pconv_ac)
    _PM.add_setpoint!(sol, pm, "convdc", "qconv", :qconv_ac)
    _PM.add_setpoint!(sol, pm, "convdc", "pdc", :pconv_dc)
    _PM.add_setpoint!(sol, pm, "convdc", "iconv", :iconv_ac)
    _PM.add_setpoint!(sol, pm, "convdc", "ptf_fr", :pconv_tf_fr)
    _PM.add_setpoint!(sol, pm, "convdc", "ptf_to", :pconv_tf_to)
    add_dcconverter_voltage_setpoint(sol, pm)
end

function add_dcgrid_flow_setpoint(sol, pm::_PM.AbstractPowerModel)
    # check the branch flows were requested
    _PM.add_setpoint!(sol, pm, "branchdc", "pf", :p_dcgrid,  status_name="status",  var_key = (idx,item) -> (idx, item["fbusdc"], item["tbusdc"]))
    _PM.add_setpoint!(sol, pm, "branchdc", "pt", :p_dcgrid;  status_name="status",  var_key = (idx,item) -> (idx, item["tbusdc"], item["fbusdc"]))
end

function add_dcbranch_losses(sol, pm::_PM.AbstractPowerModel)
    for (i, branchdc) in sol["branchdc"]
        pf = branchdc["pf"]
        pt = branchdc["pt"]
        sol["branchdc"]["$i"]["ploss"] = pf + pt
    end
end

function add_dcconverter_losses(sol, pm::_PM.AbstractPowerModel)
    for (i, convdc) in sol["convdc"]
        pf = convdc["pdc"]
        pt = convdc["pgrid"]
        ptf_fr = convdc["ptf_fr"]
        ptf_to = convdc["ptf_to"]
        sol["convdc"]["$i"]["ploss_tot"] = pf + pt
        sol["convdc"]["$i"]["ploss_tf"] = ptf_fr + ptf_to
        sol["convdc"]["$i"]["ploss_conv"] = sol["convdc"]["$i"]["ploss_tot"] - sol["convdc"]["$i"]["ploss_tf"] # TODO update later with filter etc..
    end
end

function add_dc_bus_voltage_setpoint(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "busdc", "vm", :vdcm, status_name="Vdc", inactive_status_value = 4)
end

function add_dcconverter_voltage_setpoint(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc", "vmconv", :vmc, status_name="islcc", inactive_status_value = 0)
    _PM.add_setpoint!(sol, pm, "convdc", "vaconv", :vac, status_name="islcc", inactive_status_value = 0)
    _PM.add_setpoint!(sol, pm, "convdc", "vmfilt", :vmf, status_name="islcc", inactive_status_value = 0)
    _PM.add_setpoint!(sol, pm, "convdc", "vafilt", :vaf, status_name="islcc", inactive_status_value = 0)
end

function add_dcconverter_firing_angle(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc", "phi", :phiconv, status_name="islcc", inactive_status_value = 0)
end

function get_solution_acdc_ne(pm::_PM.AbstractPowerModel, sol::Dict{String,Any})
    _PM.add_setpoint_bus_voltage!(sol, pm)
    _PM.add_setpoint_generator_power!(sol, pm)
    _PM.add_setpoint_branch_flow!(sol, pm)
    _PM.add_setpoint_dcline_flow!(sol, pm)
    if haskey(pm.data, "convdc")
        PowerModelsACDC.add_dc_bus_voltage_setpoint(sol, pm)
        PowerModelsACDC.add_dcconverter_setpoint(sol, pm)
        PowerModelsACDC.add_dcgrid_flow_setpoint(sol,pm)
        PowerModelsACDC.add_dcbranch_losses(sol, pm)
        PowerModelsACDC.add_dcconverter_losses(sol, pm)
        PowerModelsACDC.add_dcconverter_voltage_setpoint(sol, pm)
        PowerModelsACDC.add_dcconverter_firing_angle(sol, pm)
    elseif haskey(pm.data, "nw")
        if haskey(pm.data["nw"]["1"], "convdc")
            PowerModelsACDC.add_dc_bus_voltage_setpoint(sol, pm)
            PowerModelsACDC.add_dcconverter_setpoint(sol, pm)
            PowerModelsACDC.add_dcgrid_flow_setpoint(sol,pm)
            PowerModelsACDC.add_dcbranch_losses(sol, pm)
            PowerModelsACDC.add_dcconverter_losses(sol, pm)
            PowerModelsACDC.add_dcconverter_voltage_setpoint(sol, pm)
            PowerModelsACDC.add_dcconverter_firing_angle(sol, pm)
        end
    end
    if haskey(pm.data, "convdc_ne")
        add_dc_bus_voltage_setpoint_ne(sol, pm)
        add_dcconverter_setpoint_ne(sol, pm)
        add_dcgrid_flow_setpoint_ne(sol,pm)
        add_dcbranch_losses_ne(sol, pm)
        add_dcconverter_losses_ne(sol, pm)
    elseif haskey(pm.data, "nw")
        if haskey(pm.data["nw"]["1"], "convdc_ne")
            add_dc_bus_voltage_setpoint_ne(sol, pm)
            add_dcconverter_setpoint_ne(sol, pm)
            add_dcgrid_flow_setpoint_ne(sol,pm)
            add_dcbranch_losses_ne(sol, pm)
            add_dcconverter_losses_ne(sol, pm)
        end
    end
    return sol
end

function add_dc_bus_voltage_setpoint_ne(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "busdc_ne", "vm", :vdcm_ne; status_name = "Cdc", inactive_status_value = 5)
end

function add_dcconverter_setpoint_ne(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "pgrid", :pconv_tf_fr_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "qgrid", :qconv_tf_fr_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "pconv", :pconv_ac_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "qconv", :qconv_ac_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "pdc", :pconv_dc_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "iconv", :iconv_ac_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "ireactor", :irc_sq_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "iconv_sq", :iconv_ac_sq_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "ppr_fr", :pconv_pr_fr_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "ptf_to", :pconv_tf_to_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "isbuilt", :conv_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "itr", :itf_sq_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "ipr", :irc_sq_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "pfilt_to", :pconv_tf_to_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "pfilt_fr", :pconv_pr_fr_ne)

    add_dcconverter_voltage_setpoint_ne(sol, pm)
    add_dcconverter_firing_angle_ne(sol, pm)
end

function add_dcconverter_voltage_setpoint_ne(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vmconv", :vmc_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vaconv", :vac_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vmfilt", :vmf_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vafilt", :vaf_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "wf", :wf_ac_ne)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "w_du", :w_du)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "wc", :wc_ac_ne)
end

function add_dcgrid_flow_setpoint_ne(sol, pm::_PM.AbstractPowerModel)
    # check the branch flows were requested
    _PM.add_setpoint!(sol, pm, "branchdc_ne", "pf", :p_dcgrid_ne; status_name="status",  var_key = (idx,item) -> (idx, item["fbusdc"], item["tbusdc"]))
    _PM.add_setpoint!(sol, pm, "branchdc_ne", "pt", :p_dcgrid_ne; status_name="status",  var_key = (idx,item) -> (idx, item["tbusdc"], item["fbusdc"]))
    _PM.add_setpoint!(sol, pm, "branchdc_ne", "isbuilt", :branchdc_ne)
    _PM.add_setpoint!(sol, pm, "branchdc_ne", "wdc_to", :wdc_du_to; status_name="status")
end

function add_dcbranch_losses_ne(sol, pm::_PM.AbstractPowerModel)
    for (i, branchdc) in sol["branchdc_ne"]
        pf = branchdc["pf"]
        pt = branchdc["pt"]
        sol["branchdc_ne"]["$i"]["ploss"] = pf + pt
    end
end

function add_dcconverter_losses_ne(sol, pm::_PM.AbstractPowerModel)
    for (i, convdc) in sol["convdc_ne"]
        pdc = convdc["pdc"]
        pgrid = convdc["pgrid"]
        ptf_to = convdc["ptf_to"]
        pconv = convdc["pconv"]
        sol["convdc_ne"]["$i"]["ploss_tot"] = pdc - pgrid
        sol["convdc_ne"]["$i"]["ploss_tf"] = pgrid + ptf_to
        sol["convdc_ne"]["$i"]["ploss_re"] = - ptf_to + pconv # ptf_to + ppr_fr = 0
        sol["convdc_ne"]["$i"]["ploss_conv"] =  pconv + pdc
    end
end

function add_dcconverter_firing_angle_ne(sol, pm::_PM.AbstractPowerModel)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "phi", :phiconv_ne)
end


text:
module TimescaleDBs

using HTTP, JSON2
using Dates
using DataFrames

using ..AbsTypes
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!

export fit, fit!,transform, transform!
export TimescaleDB

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.
mutable struct TimescaleDB <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function TimescaleDB(args=Dict())
      default_args = Dict(
         :uri => "http://localhost:3000",
         :db => "dateval",
         :dbusername =>"postgres",
         :dbpassword => "mysecretpassword",
         :query => "select=date,value"
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

function fit!(tdb::TimescaleDB, features::DataFrame=DataFrame(), labels::Vector=[])::Nothing
   (features == DataFrame() && labels == [])  || throw(ArgumentError("features and labels should be empty because data are from http request"))
   uri = tdb.model[:uri]; db = tdb.model[:db]; query = tdb.model[:query]
   (uri != "" && query != "" && db != "") || error("missing uri/query/db")
   return nothing
end

function fit(tdb::TimescaleDB, features::DataFrame=DataFrame(), labels::Vector=[])::TimescaleDB
   fit!(tdb,features,labels)
   return deepcopy(tdb)
end

function transform!(tdb::TimescaleDB, features::DataFrame=DataFrame())::DataFrame
   features == DataFrame()  || throw(ArgumentError("features should be empty because data are from http request"))
   uri = tdb.model[:uri]; db = tdb.model[:db]; query = tdb.model[:query]
   (uri != "" && query != "" && db != "") || error("missing uri/query/db")
   payload = uri*"/"*db*"?"*query
   req = HTTP.get(payload)
   req.status == 200 || error(req.status)
   body = String(req.body)
   df = body |> JSON2.read |> DataFrame
   ncol(df) == 2 || error("data should have 2 columns")
   nrow(df) > 0 || return DataFrame()
   rename!(df,names(df)[1] => :Date, names(df)[2] => :Value)
   df.Date = DateTime.(df.Date)
   eltype(df.Value) <: Number || throw(ArgumentError("values in second column are not numeric"))
   return df
end

function transform(tdb::TimescaleDB, features::DataFrame=DataFrame())::DataFrame
   return transform!(tdb,features)
end

end


text:
"""
    GMC_NS_Model_Record

Abstract supertype for records of serialised 'GMC_NS_Model's. Subtyped record structs must implement the following fields as their interface with 'GMC_NS':

    struct Minimal_Record
        trajectory::Integer #id assigned at construction, identifies trajectory
        i::Integer #id integer assigned at construction, identifies position on trajectory
        pos::Vector{<:AbstractFloat} #hypersphere position of model-particle
        path::String #model is serialized to ensembledir/trajectory.i
        log_Li::AbstractFloat #log likelihood calculated by 'GMC_NS_Model' constructor
    end
"""
abstract type GMC_NS_Model_Record end

"""
    GMC_NS_Model_Record

Abstract supertype for models in 'GMC_NS_Ensemble's. Subtyped model structs must NOT be mutable, and must not have inner constructors that make field-by-field construction inaccessible. Subtyped model structs must implement the following fields as their interface with 'GMC_NS', in order, and all additional fields must follow:

    mutable struct Minimal_Model
        trajectory::Integer #id assigned at construction, identifies trajectory
        i::Integer #id integer assigned at construction, identifies position on trajectory

        θ::Vector{<:AbstractFloat} #model's parameter Vector
        log_Li<:AbstractFloat #model's likelihood, calculated by constructor

        pos::Vector{<:AbstractFloat}
        v::Vector{<:AbstractFloat} #model's velocity in hypersphere coords
    end
"""
abstract type GMC_NS_Model end