Using type parameters in the types of the fields

General Usage
,
1

I’m not a data structures person, so the right answer here might be “you’re doing it wrong, try this instead”; but I’ll start with the immediate question I have:

in a recursive type, is it not possible to have the type parameters of the fields be related to the type parameters of the struct itself? for example, if one were to try to define:

struct Simplex{T<:Any, dim, size}
    x::Union{NTuple{size, Simplex{T, dim-1}}, T}
end

you get no method matching -(::Typevar, ::Int64).

what I'm actually trying to do

I have data which is ‘shaped’ like a n-dimensional regular simplex, by which I mean I often want the data which is along some specific edge, face, volume, hyper-volume, etc of the simplex. eg, say the data is 1:10 and the dimension is 2 (so size is 4), I would want to represent it like

(
    (1),
    (2, 3),
    (4, 5, 6),
    (7, 8, 9, 10)
)

since I’ll be needing last() (7, 8, 9, 10), last.() (1, 3, 6, 10), and first.() (1, 2, 4, 7) a lot.
or with the same data if the dimension (and size) was 3, I would want it to look like

(
    ((1)),
    ((2), (3, 4)),
    ((5), (6, 7), (8, 9, 10))
)

and the pieces I’d need are last() ((5), (6,7), (8,9,10)), last.() ((1), (3,4), (8,9,10)), last..() ((1), (2,4), (5,7,10)), and first..() ((1), (2,3), (5,6,8)). (if you’re really good at visualizing this kind of thing, you’ll notice those are the four faces on this 3d triangular pyramid)

most of them are 4- or 5- d, where faking it with an array starts getting really confusing; but they can be arbitrarily big (seen up to 10-d so far)

I tried defining it as

b(size, dim) = prod(size:dim+size-1)/prod(1:dim)

struct Simplex{T<:Any, dim, size}
    x::Union{NTuple{size, Simplex{T, dimMinusOne}}, T} where dimMinusOne
    Simplex{T, dim, size}(x::Vector{T}) where {T} = dim==1 ?
        new{T, 1, size}(ntuple(n->Simplex{T}(x[n]), size)) :
        new{T, dim, size}(ntuple(n->Simplex{T, dim-1, n}(x[b(n-1,dim)+1:b(n,dim)]), size))
    Simplex{T}(x::T) where T = new{T, 0, 1}(x)
end

but Simplex{Int, 2, 2}([1, 2, 3]) says there’s no method.

2

It’s not possible to do math on typevars (though just today I see Proof of concept: Allow some limited arithmetic on type vars - Pull Request #62707 - JuliaLang/julia - GitHub). This is why StaticArray’s SMatrix is parameterized like SMatrix{M, N, T, MtimesN} with a “redundant” fourth parameter that is simply required to be the product M*N (the length, which parameterizes the NTuple used for storage).

In your case you can do something similar:

struct Simplex{T<:Any, dim, size, dim1}
    x::Union{NTuple{size, Simplex{T, dim1}}, T}
end
# use inner constructor to set (and require) dim1 = max(dim-1,0)

Do note that your inner simplex is not fully typed (missing size and its own dim1) so you may have a bit of boxing and type instability here.

Personally, I would use a design more like

struct Simplex{T, dim, size, SubSimplex}
    x::NTuple{size, SubSimplex}
end

where SubSimplex is required to be either <:Simplex{T, dim-1} or T. You could use the Union{..., T} like before, but instead I would simply set size=1, SubSimplex=T at the bottom level.

Although I would probably go further and not have a naked T at all, simply define the (0 or 1, I haven’t thought about it) dimensional Simplex{T} to follow all the relevant rules of a scalar.

3

yeah, it’s the 0-d simplex that’s a point; 1-d is a line (vector). first two dimensions are exactly the same as arrays. all three possible n-d regular polytopes are the same in 0- and 1- d.

this sounds like exactly what I want, but still getting a a syntax error trying to define it:

struct Simplex{T, dim, size, SubSimplex} where {SubSimplex <: Simplex{T, dim-1}}
    x::NTuple{size, SubSimplex}
end
ERROR: syntax: invalid type signature


struct Simplex{T, dim, size, SubSimplex} where {SubSimplex isa Simplex{T, dim-1}}
    x::NTuple{size, SubSimplex}
end
ERROR: syntax: invalid type signature

the reason I tried to put a union with the naked T is because if the NTuple is made of Simplex and Simplex is made of NTuple then I don’t know how to stop the infinite loop; if there’s a better way around that, I’m all ears ^^

4

Just do the definition as I wrote it. You can’t enforce the where in the struct definition. The place to enforce this is in an inner constructor.

Here’s a sketch (that almost certainly has multiple errors – I didn’t even check whether this was syntactically correct):

struct Simplex{T, dim, size, SubSimplex}
    x::NTuple{size, SubSimplex}
    function Simplex{T, dim, size, SubSimplex}(subsimplices::NTuple{size, SubSimplex}) where {T, dim, size, SubSimplex}
        SubSimplex <: Simplex{T, dim-1} || error("subsimplices of Simplex{T, dim} must be Simplex{T, dim-1}") # don't allow construction of nonsense simplex
        return new{T, dim, size, SubSimplex}(subsimplices)
    end
end

dimension(::Type{<:Simplex{<:Any, dim}}) where dim = dim
dimension(x::Simplex) = dimension(typeof(x))
Base.eltype(::Type{<:Simplex{T}}) where T = T
Base.eltype(x::Simplex) = Base.eltype(typeof(x))

function Simplex(subsimplices::NTuple{size, SubSimplex}) where {size, SubSimplex}
    # fill parameters to call inner constructor
    return Simplex{eltype(SubSimplex), dimension(SubSimplex)+1, size, SubSimplex}(subsimplices)
end

For a better example, maybe look at how StaticArrays does it.

5

Related links:

6

It’s not that it’s not possible to do math on typevars, exactly, because we can enforce invariants in constructors. It’s that we can dynamically redefine that math in the general case. If dim-1 can be different values across different world ages for the same dim value, then dim-1 must be an independent parameter to distinguish the different concrete subtypes of the same parametric type, no matter how redundant it often appears. I interpreted freezing the world age to be part of the follow-up draft (types: Support computed field types via world-age-detached field gens by Keno · Pull Request #62715 · JuliaLang/julia · GitHub), though it’s not clear to me how that would stay mentally consistent across technically distinct but related parametric types bound in different world ages, including redefined structs. Then again, since the requests for type parameter computation are usually for basic arithmetic on a fixed set of concrete types, maybe parameter arithmetic being extended across invisible world ages is a more theoretical than practical problem.

7

thanks, I see I was trying to put too much into the struct definition. hadn’t thought of letting the struct itself be more permissive and put the restrictions in a separate constructor function. still wasn’t able to eliminate the Union{} from my definition; how hard should I try to find a way to take it out?

struct Simplex{dim, eltype, size}
    x::Union{NTuple{size, Simplex}, eltype}
end

b(size, dim) = prod(size:dim+size-1)÷prod(1:dim)
function Simplex{dim}(vals::Vector{eltype}) where {dim, eltype}
    if dim==0 
        @assert length(vals)==1 "0-d simplicies must have exactly one element"
        Simplex{0, eltype, 1}(vals[])
    else
        len, elms = 0, 0
        while elms<length(vals)
            len+=1
            elms=b(len, dim)
        end
        @assert b(len, dim) == length(vals) "wrong number of vals. expected $(b(len, dim)), but only got $(length(vals))"
        Simplex{dim, eltype, len}(Simplex{dim-1}.(ntuple(n->vals[b(n-1,dim)+1:b(n,dim)], len)))
    end
end
now I'm off to figure out how to make distribution put the result back into my type instead of eating away at my dimensions... 
dim(i::Simplex) = typeof(i).parameters[1]
leaftype(i::Simplex) = typeof(i).parameters[2]
Base.summary(io::IO, x::Simplex) = sprint(
    print,
    !iszero(dim(x)) ? "$(length(x))-element $(dim(x))" : '0', "-simplex with eltype ", leaftype(x);
    context=io
)
import REPL: show_repl
function show_repl(io::IO, m::MIME"text/plain", x::Simplex)
    println(io, summary(io, x), ':')
    iszero(dim(x)) ? (print(io, ' '); show(io, m, x)) :
    for i in x
        print(io, ' ')
        show(io, m, i)
        i === x[end] || println(io)
    end
end
function Base.show(io::IO, m::MIME"text/plain", x::Simplex)
    if dim(x) == 0
        show(io, m, x.x)
    else
        print(io, '(')
        for i in x
            show(io, m, i)
            i === x[end] || print(io, ", ")
        end
        print(io, ')')
    end
end

Base.length(iter::Simplex) = length(iter.x)
Base.getindex(iter::Simplex, i::Int) = getindex(iter.x, i)
Base.firstindex(iter::Simplex) = firstindex(iter.x)
Base.lastindex(iter::Simplex) = lastindex(iter.x)
Base.iterate(iter::Simplex) = iter.x[1], 1
Base.iterate(iter::Simplex, state) = length(iter)>state ? (iter.x[state+1], state+1) : nothing

julia> a = Simplex{4}([1:15;])
3-element 4-simplex with eltype Int64:
 (((1)))
 (((2)), ((3), (4, 5)))
 (((6)), ((7), (8, 9)), ((10), (11, 12), (13, 14, 15)))

julia> last(a)
3-element 3-simplex with eltype Int64:
 ((6))
 ((7), (8, 9))
 ((10), (11, 12), (13, 14, 15))

julia> last.(a)
3-element Vector{Simplex{2, Int64}}:
 ((1))
 ((3), (4, 5))
 ((10), (11, 12), (13, 14, 15))

julia> last'²(a)
3-element Vector{Vector}:
 Simplex{1, Int64, 1}[(1)]
 Simplex{1, Int64}[(2), (4, 5)]
 Simplex{1, Int64}[(6), (8, 9), (13, 14, 15)]

julia> last'³(a)
3-element Vector{Vector{Vector{Simplex{0, Int64, 1}}}}:
 [[1]]
 [[2], [3, 5]]
 [[6], [7, 9], [10, 12, 15]]

julia> first'³(a)
3-element Vector{Vector{Vector{Simplex{0, Int64, 1}}}}:
 [[1]]
 [[2], [3, 4]]
 [[6], [7, 8], [10, 11, 13]]
8

Yeah looks like I hadn’t thought things all the way through when I suggested you could eliminate the Union. An obvious way to remove it isn’t coming to mind. The best I have is something like the fully-parameterized types I was proposing above:

struct Simplex{dim, eltype, SubsimplicesType}
    x::SubsimplicesType
end

and now Subsimplices can be either the NTuple of Simplex (now able to include all its type parameters, though your types get increasingly complicated with dimension) or the bottom eltype. With this, you could eliminate the eltype parameter too but perhaps you want that for dispatch or easy access.

9

would doing that improve type stability or something? I really don’t have a good mental model for how different type definitions effect performance, but that sounds to me like it’d just be renaming Union{...} to SubSimpliciesType, which I would think would have the same performance?

the solution to broadcasting was way too simple for it to have taken 6 hours...

why does our documentation not make a single mention of the one actually useful broadcasting function?

Broadcast.broadcastable(s::Simplex) = s
Broadcast.BroadcastStyle(::Type{<:Simplex}) = Broadcast.Style{Simplex}()
function Broadcast.materialize(B::Broadcast.Broadcasted{Broadcast.Style{Simplex}})
    B = Broadcast.flatten(B)
    args = ntuple(n->B.args[n] isa Simplex ? (i for i=B.args[n]) : B.args[n], length(B.args))
    vec = B.f.(args...)
    Simplex{dim(vec[1])+1, leaftype(vec[1]), length(vec)}(ntuple(n->vec[n],length(vec)))
end

julia> a = Simplex{4}([1:15;])
3-element 4-simplex with eltype Int64:
 (((1)))
 (((2)), ((3), (4, 5)))
 (((6)), ((7), (8, 9)), ((10), (11, 12), (13, 14, 15)))

julia> last(a)
3-element 3-simplex with eltype Int64:
 ((6))
 ((7), (8, 9))
 ((10), (11, 12), (13, 14, 15))

julia> last.(a)
3-element 3-simplex with eltype Int64:
 ((1))
 ((3), (4, 5))
 ((10), (11, 12), (13, 14, 15))

julia> last'²(a)
3-element 3-simplex with eltype Int64:
 ((1))
 ((2), (4, 5))
 ((6), (8, 9), (13, 14, 15))

julia> last'³(a)
3-element 3-simplex with eltype Int64:
 ((1))
 ((2), (3, 5))
 ((6), (7, 9), (10, 12, 15))

julia> first'³(a)
3-element 3-simplex with eltype Int64:
 ((1))
 ((2), (3, 4))
 ((6), (7, 8), (10, 11, 13))
10

The issue with NTuple{size, Simplex} is that Simplex is basically Simplex{<:Any, <:Any, <:Any}, which is to say that the compiler knows virtually nothing about what its type parameters or field types are. You might get a partial improvement from instead defining NTuple{size, Simplex{Any, eltype, Any} so that it knows one of the parameters.

In other words, any access is to a value of unknowable-at-compile-time type, which means type instability and dynamic dispatch which has some fixed overhead (which can add up to a lot or a little performance overhead, depending on your use case). It also means boxing of variables. The compiler not knowing size probably makes your NTuple{size, X} no better than Vector{X} (probably worse, were I to guess).

Using a type parameter lets the parameters be concrete, which resolves this issue. But there exist situations where it won’t really be possible to make these concrete regardless (in which case the parameters will be non-concrete and won’t really help).

11

I had an earlier version that looks like this:

struct Simplex{type, dim, eltype, size}
    x::Union{NTuple{size, <:type}, eltype}
end

b(size, dim) = prod(size:dim+size-1)÷prod(1:dim)
function Simplex{dim}(vals::Vector{eltype}) where {dim, eltype}
    let len=0, elms=0
        while elms<length(vals)
        len+=1
        elms=b(len, dim)
    end
    b(len, dim) == length(vals) || error("wrong number of vals. expected $(b(len, dim)), but got $(length(vals))")
    dim==0 ?
        Simplex{Simplex{eltype}, 0, eltype, 1}(vals[]) :
        Simplex{Simplex{dim==1 ? eltype : Simplex, dim-1, eltype}, dim, eltype, len}(Simplex{dim-1}.(ntuple(n->vals[b(n-1,dim)+1:b(n,dim)], len)))
    end
end

which works for 0- and 1- d, but since <: doesn’t propagate it breaks with higher dimensions. (since Tuple{Int} <: Tuple{Real} == true I’m sure there’s a way to fix that but I gave up.) does that help anything? it has a type that would look like

julia> typeof(Simplex{0}([1]))
Simplex{Int64, 0, Int64, 1}

julia> typeof(Simplex{1}([1:10;]))
Simplex{Simplex{Int64, 0, Int64}, 1, Int64, 10}

julia> typeof(Simplex{2}([1:10;]))
Simplex{Simplex{Simplex, 1, Int64}, 2, Int64, 4}

julia> typeof(Simplex{3}([1:10;]))
Simplex{Simplex{Simplex, 2, Int64}, 3, Int64, 3}

julia> typeof(Simplex{4}([1:15;]))
Simplex{Simplex{Simplex, 3, Int64}, 4, Int64, 3}

size would be especially hard to inform the compiler about, since it’s different per each inner element; I’d guess you can at the most tell it the number of elements in the outer layer.

wouldn’t

struct Simplex{dim, eltype, SubsimplicesType}
    x::SubsimplicesType
end

just require a

struct SubsimpliciesType{dim, eltype, size}
    x::Union{NTuple{size, Simplex}, eltype}
end

itself? how does it not just push the problem one alias down?

12

since <: doesn’t propagate

Note that using <:X in a struct field is probably seldom-better than Any (maybe unless X is concrete, but even then it might not be that smart).

The idea is that Simplex{dim, eltype, size, SubsimplicesType} would have

SubsimplicesType == Simplex{dim-1, eltype, size1,
  Simplex{dim-2, eltype, size2,
    Simplex{dim-3, eltype, size3,
      ...
        Simplex{0, eltyle, sizeN, eltype}
      ...
    }}}

Note that if you can start with a concrete base element and construct the data contents of the Simplex, the compiler should be able to construct these types for you.

size would be especially hard to inform the compiler about, since it’s different per each inner element

Does the above address your size issue, allowing each level to be different? Or can each element of the NTuple have a different size?

If the case is different-per-element, it’s not really a compile-time-knowable value and NTuple is a pretty bad container for it. A Vector{T} will outperform a NTuple{Any, T} (i.e., with uninferred length) in almost all contexts.

If it’s per-level, I think the above will work for you (up to a moderate number of dimensions, anyway – at larger sizes a Vector might still outperform NTuple).

Note that in either case the complexity of the type is dim types “deep.” This will get ugly and can eventually become a problem in high dimensions (is probably okay for at least 10-30 dimensions).

I haven’t worked with simplices in long enough that I really don’t remember how they operate. From your examples, it looks like you can easily “flatten” all the data for an arbitrary-dimension simplex. In which case an alternative (and probably better) storage format might be something like

struct Simplex{dim, elT}
  data::Vector{elT}
end
#=
# alternative: parameterize the container type as well
# useful if subsimplex data has contiguous indices so `SubArray` (aka `view`) is useful
struct Simplex{dim, elT, D}
  data::D
end
=#

# could define this as Base.getindex if you want simplex[i] syntax
function getsubsimplex(s::Simplex{dim}, index) where {dim}
  # construct and return the `Simplex{dim-1, ...}` (or scalar) at `index`
end

where you simply construct each subsimplex on-demand from its parent.

13

each element of the NTuple is a strictly different size (except for, 0- and 1- d, since 0-d always has size of 1). using the examples I had put in the original post (looks like I put the examples under the [details] fold. oops.):

2-d (triangle):

sim2 = ( #  <- top-level Simplex, currently of type `Simplex{2, Int, 4}`
    (1),    #  <- first element, type `Simplex{1, Int, 1}`
    (2, 3),    #  <- second element, `Simplex{1, Int, 2}`
    (4, 5, 6),    #  <- `Simplex{1, Int, 3}`
    (7, 8, 9, 10)    #  <- `Simplex{1, Int, 4}`
)

where the faces are last(sim2) == (7, 8, 9, 10), last.(sim2) == (1, 3, 6, 10), and first.(sim2) == (1, 2, 4, 7).

and to give a feel for how it generalizes, 3-d (a tetrahedron):

tetra = ( # has type `Simplex{3, Int, 3}`
    ((1)),    # `Simplex{2, Int, 1}`
    ((2), (3, 4)),    # `Simplex{2, Int, 2}`
    ((5), (6, 7), (8, 9, 10))    # `Simplex{2, Int, 3}`
)

with triangular faces:
last(tetra) == ((5), (6,7), (8,9,10)) (corners 5, 8, 10),
last.(tetra) == ((1), (3,4), (8,9,10)) (corners 1, 8, 10),
last..(tetra) == ((1), (2,4), (5,7,10)) (corners 1, 5, 10), and
first..(tetra) == ((1), (2,3), (5,6,8)) (corners 1, 5, 8)

in other words, the first element is always of size 1, and the nth element is of size n.

I have something similar to Base.getindex(s, i...) = [s=s.x[k] for k=i][end] for indexing currently; so tetra[3, 2, 2] is 7.

yeah, if you know the dimension you could enumerate the indices and store it as a vector. I mean, my initializer takes in {dim} and vals::Vector{eltype} so if I wasn’t able to figure out the shape from that info I’d be in trouble.
was kinda hoping to frontload computation in initialization so that getindex is easy; but it would be good to actually test if that does help me. it’s also slightly complicated by the fact that I’m doing a bit of data reorganization in my actual Simplex{dim}() which is basically free when defined this way, will have to think on it longer to see if I can make getindex do that same reorganization on a vector for ‘free’ as well.

14

ok, after spending way too long figuring out how correctly implement broadcasting for the single-vector storage method, I can now compare the three methods:

the one I’ve been using above, x::Union{NTuple{size, Simplex}, eltype}}

julia> @benchmark n_faces($2,$a)
BenchmarkTools.Trial: 29 samples with 1 evaluation per sample.
 Range (min … max):  152.015 ms … 240.883 ms  ┊ GC (min … max): 0.00% … 12.40%
 Time  (median):     166.473 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   174.199 ms ±  21.629 ms  ┊ GC (mean ± σ):  2.42% ±  4.36%

   █   ▃  ▃█ █        ▃    ▃                ▃
  ▇█▇▇▇█▇▁██▁█▁▁▇▁▇▁▁▁█▁▁▁▁█▁▇▁▁▁▁▁▁▁▇▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
  152 ms           Histogram: frequency by time          241 ms <

 Memory estimate: 11.84 MiB, allocs estimate: 279768.

naive tuple → vec replacement (x::Union{Vector{Simplex}, eltype})

julia> @benchmark n_faces($2,$a)
BenchmarkTools.Trial: 33 samples with 1 evaluation per sample.
 Range (min … max):  139.688 ms … 189.200 ms  ┊ GC (min … max): 0.00% … 17.17%
 Time  (median):     147.847 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   155.865 ms ±  16.149 ms  ┊ GC (mean ± σ):  4.47% ±  7.39%

       ▂  █▂▅                                         ▂
  █▅▁▁▁█▅▅███▁▅▁▅▁▁▁▁▅▁▅▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▅▅▁▁▁▁█ ▁
  140 ms           Histogram: frequency by time          189 ms <

 Memory estimate: 17.38 MiB, allocs estimate: 447322.

stored as a single vector (x::Vector{eltype})

julia> @benchmark n_faces($2,$a)
BenchmarkTools.Trial: 25 samples with 1 evaluation per sample.
 Range (min … max):  180.130 ms … 264.397 ms  ┊ GC (min … max): 0.00% … 16.44%
 Time  (median):     189.923 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   204.319 ms ±  26.661 ms  ┊ GC (mean ± σ):  6.05% ±  8.31%

     ▁█▄▁▁  ▁                            ▄
  ▆▁▁█████▆▁█▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▆▁▁▁▆▁▁▁▁▁▆▁▁▁▁▁▆ ▁
  180 ms           Histogram: frequency by time          264 ms <

 Memory estimate: 35.84 MiB, allocs estimate: 803199.

I think most of the minus from storing it as a single vector is because the n_faces command needs a frankly insane amount of broadcasting, and if it’s stored as a single vector that means we have to keep splitting up and re-collecting it a bunch.

as expected, it makes creation a flash:

julia> @benchmark a = Simplex{10}($[1:3003;]) # tuples
BenchmarkTools.Trial: 341 samples with 1 evaluation per sample.
 Range (min … max):  12.297 ms … 75.592 ms  ┊ GC (min … max): 0.00% … 75.56%
 Time  (median):     12.861 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   14.671 ms ±  7.724 ms  ┊ GC (mean ± σ):  6.99% ± 10.91%

  █▄
  ███▄▅▄▅▇█▅▅▁▅▆▄▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▄▄ ▆
  12.3 ms      Histogram: log(frequency) by time      62.8 ms <

 Memory estimate: 2.86 MiB, allocs estimate: 94054.

julia> @benchmark a = Simplex{10}($[1:3003;]) # vector of simplex
BenchmarkTools.Trial: 300 samples with 1 evaluation per sample.
 Range (min … max):  13.800 ms … 66.840 ms  ┊ GC (min … max):  0.00% … 76.27%
 Time  (median):     14.609 ms              ┊ GC (median):     0.00%
 Time  (mean ± σ):   16.685 ms ±  9.382 ms  ┊ GC (mean ± σ):  10.80% ± 14.38%

  █▆▁▂
  ████▇▄▁▁▅▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▁▆▁▆ ▆
  13.8 ms      Histogram: log(frequency) by time      65.6 ms <

 Memory estimate: 4.86 MiB, allocs estimate: 134381.

julia> @benchmark a = Simplex{10}($[1:3003;]) # single vector
BenchmarkTools.Trial: 10000 samples with 216 evaluations per sample.
 Range (min … max):  344.713 ns …  2.054 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     349.875 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   355.530 ns ± 32.164 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▇█▇▅▄▂▁ ▁▂▂▂▂▁                                              ▂
  █████████████████▇▇▆▆▆▆▆▆▆▅▅▄▄▆▅▆▆▅▅▄▅▃▆▅▄▅▅▄▅▄▄▄▄▂▃▄▄▃▄▄▂▂▄ █
  345 ns        Histogram: log(frequency) by time       457 ns <

 Memory estimate: 16 bytes, allocs estimate: 1.

and makes accessing elements slower:

julia> @benchmark $a[$5,$5,$5,$3,$2,$2] # tuple
BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample.
 Range (min … max):  3.570 μs …  13.700 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.721 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.855 μs ± 504.136 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▄▇█▇▅▂▂▃▂▂                                                 ▂
  ▆███████████▇▇▇▇▇▇█▇▇▇▇▇▇▇▇▇▇▇▆▆▅▅▆▆▆▆▆▅▆▆▆▅▆▅▄▇▆▆▇▆▇▆▆▆▆▅▅ █
  3.57 μs      Histogram: log(frequency) by time      5.93 μs <

 Memory estimate: 128 bytes, allocs estimate: 2.

julia> @benchmark $a[$5,$5,$5,$3,$2,$2] # vector of simplex
BenchmarkTools.Trial: 10000 samples with 26 evaluations per sample.
 Range (min … max):  930.385 ns … 126.352 μs  ┊ GC (min … max): 0.00% … 96.69%
 Time  (median):     961.385 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.014 μs ±   1.264 μs  ┊ GC (mean ± σ):  1.20% ±  0.97%

  ▂▇█▇▆▅▄▃▃▂▁▁  ▁▁▁▁                                     ▁▂▂▁   ▂
  █████████████████████▇▇▇▅▆▆▆▆▆▅▇▆▆▆▆▅▄▆▆▅▄▅▃▅▅▄▄▄▁▄▃▃▃▇█████▆ █
  930 ns        Histogram: log(frequency) by time       1.54 μs <

 Memory estimate: 128 bytes, allocs estimate: 2.

julia> @benchmark $a[$5,$5,$5,$3,$2,$2] # single vector
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
 Range (min … max):  5.488 μs … 659.029 μs  ┊ GC (min … max):  0.00% … 90.60%
 Time  (median):     6.287 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   8.711 μs ±  25.734 μs  ┊ GC (mean ± σ):  19.72% ±  6.83%

     ▆██▄▁
  ▁▃▇██████▇▆▅▄▃▂▂▂▂▂▂▂▁▁▁▂▂▃▃▃▃▂▂▂▂▂▁▁▁▁▁▁▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  5.49 μs         Histogram: frequency by time          12 μs <

 Memory estimate: 18.97 KiB, allocs estimate: 33.

since I’m creating it once at the beginning and then basically working with that for the rest of the time, looks like I’m going with Union{Vector{Simplex},eltype}.

the `n_faces` command: 
using Iterators: product, repeated, reverse

rebrod(f, n) = foldl(|>, repeated(Base.BroadcastFunction, n); init=f)

"return a vector containing all n-dimensional faces of an arbitrary simplex"
function n_faces(n::Int, s::Simplex{dim, eltype, len}) where {dim, eltype, len}
    0≤n≤dim || error("n must be between 0 and dim")
    n = dim-n
    n==0 && return [s]
    ind = b(dim-n+2,n)+1
    out = Vector{Simplex{dim-n, eltype, dim==n ? 1 : len}}(undef, ind-1)
    n_face = Vector{Simplex}(undef, n-1)
    n==1 && (outvec!(out, ind, 1, dim, s), return out)
    for k in product(repeated(1:dim-n+2,n-1)...) .|> reverse
        issorted(k) || continue
        rep = findfirst(==(k[end]), k)
        i = rep==1 ? s : n_face[rep-1]
        n_face[rep:end] = [i = rebrod(last, k[j]-1)(i) for j=rep:n-1]
        ind = outvec!(out, ind, k[end], dim-n+1, i)
    end
    out
end

function outvec!(out, ind, k, l, elm)
    for i in k:l
        out[(ind-=1)] = rebrod(last, i-1)(elm)
    end
    out[(ind-=1)] = rebrod(first, l-1)(elm)
    ind
end
15

If you’re after performance, note that most of those benchmarks show a rather large number of allocations. While a few appear unavoidable, it looks like a lot are the result of type instability that could be improved.

I took a crack at a flattened storage version. Try this (and check that I got it right in higher dimensions):

struct Simplex{dim, T}
    data::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true} # view(::Vector{T}, ::UnitRange{Int64})
end

Simplex{dim}(s::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}) where {dim, T<:Real} = Simplex{dim, T}(s) # fill type parameter
Simplex{dim}(s::Vector{T}) where {dim, T<:Real} = Simplex{dim, T}(@view s[begin:end]) # convert to subarray
Simplex{dim}(x::AbstractVector{T}) where {dim, T<:Real} = Simplex{dim}(convert(Vector{T}, x)) # convert to Vector and pass to Vector method

Base.getindex(s::Simplex{0}, i::Int) = s.data[i] # TODO: should only allow i==1
function Base.getindex(s::Simplex{dim}, i::Int) where dim
    common = prod(i : i-2+dim)
    start = common * (i-1) ÷ factorial(dim) + 1 # prod(i-1 : i-2+dim) ÷ factorial(dim) + 1
    stop = common * (i-1+dim) ÷ factorial(dim) # prod(i : i-1+dim) ÷ factorial(dim)
    # TODO: better error message for out-of-range i
    return Simplex{dim-1}(view(s.data, start:stop))
end

You could further parameterize the storage type, but the SubArray-of-Vector seemed to cover your needs here.

julia> using BenchmarkTools

julia> dat10 = Vector(1:10);

julia> dat3003 = Vector(1:3003);

julia> s2 = Simplex{2}(dat10)
Simplex{2, Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

julia> s3 = Simplex{3}(dat10)
Simplex{3, Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

julia> [s2[i] for i in 1:4]
4-element Vector{Simplex{1, Int64}}:
 Simplex{1, Int64}([1])
 Simplex{1, Int64}([2, 3])
 Simplex{1, Int64}([4, 5, 6])
 Simplex{1, Int64}([7, 8, 9, 10])

julia> [s3[i] for i in 1:3]
3-element Vector{Simplex{2, Int64}}:
 Simplex{2, Int64}([1])
 Simplex{2, Int64}([2, 3, 4])
 Simplex{2, Int64}([5, 6, 7, 8, 9, 10])

julia> @benchmark Simplex{10}($dat3003)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  1.715 ns … 88.216 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.758 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.809 ns ±  1.262 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

      █
  ▄▃▂▇█▄▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▁▂▂▂▂▁▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂ ▂
  1.72 ns        Histogram: frequency by time        2.32 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> a = Simplex{10}(dat3003);

julia> @benchmark $a[$5][$5][$5][$3][$2][$2]
BenchmarkTools.Trial: 10000 samples with 995 evaluations per sample.
 Range (min … max):  26.067 ns … 356.634 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     27.665 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   30.194 ns ±  11.121 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▆▅▅▇▄▃▁▁ ▁▂▁ ▁▁                                             ▁
  ████████████████▇▇▆▆▆▅▇▇▇▆▄▆▄▇▅▅▅▅▂▅▅▅▅▅▃▅▄▄▅▆▄▄██▅▂▂▅▃▄▄▆▅▄ █
  26.1 ns       Histogram: log(frequency) by time        70 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

To get full performance on n_faces, you can hopefully adjust it to use fully-typed Vector{Simplex{dim, T, <any_other_params>}{ instead of an incompletely-typed Vector{Simplex}.

16

oh yeah, forgot about @views. slicing and dicing the indices rather than manipulating the actual memory values sounds like it certainly would be a lot more efficient. I’ll go on another quest converting the one-vector implementation I’ve got into a views version ^^

I don’t think it should be quite as bad as the recursive->vector transition since I’ve got the pattern figured out now, but may have to make a specialized version of first and last to get the broadcasting working right. we’ll see what happens.

I don’t particularly need performance that bad, but it’s a fun problem to chew on. if I didn’t at least give a good attempt trying it out, I wouldn’t be satisfied. this is ultimately something I’m doing because it’s fun, no reason to not make it performant. the benefit of personal projects is that if you fall into the trap of premature optimization you can decide after the fact that the project was actually just about making that one function really performant ^^

17

@views ended up not being helpful:

julia> @benchmark n_faces($2,$a)
BenchmarkTools.Trial: 12 samples with 1 evaluation per sample.
 Range (min … max):  354.784 ms … 741.294 ms  ┊ GC (min … max):  7.32% … 33.07%
 Time  (median):     379.121 ms               ┊ GC (median):     7.05%
 Time  (mean ± σ):   419.450 ms ± 106.107 ms  ┊ GC (mean ± σ):  10.29% ±  7.96%

    █▃
  ▇▇██▁▁▇▁▁▇▁▇▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
  355 ms           Histogram: frequency by time          741 ms <

 Memory estimate: 100.07 MiB, allocs estimate: 1891567.

(with type):

struct Simplex{dim, eltype, size, rev} # rev is a bool. needed for reasons...
    x::SubArray{eltype, 1, Vector{eltype}, Tuple{Vector{Int}}, false}
end

an arbitrary face is not necessarily contiguous in memory (eg, the left side of a triangle, [1, 2, 4, 7]); so using views turns out to be basically the same as working with the data itself, except you’re splitting and re-joining the indices vector instead of the data vector. storing twice as many vectors where one is basically a vector of pointers into the other doesn’t improve anything, it just has to follow an extra pointer every time it tries to grab a value now