Submanifold index arrays should be typed

[eschnett]

trafficstars

I received this error:

Linear operators for D=1: Error During Test at /Users/eschnett/src/jl/RayTraceGR/test/runtests.jl:52
  Got exception outside of a @test
  MethodError: no method matching one(::Type{Any})
  Closest candidates are:
    one(::Type{Union{Missing, T}}) where T at missing.jl:105
    one(!Matched::Type{Missing}) at missing.jl:103
    one(!Matched::BitArray{2}) at bitarray.jl:422
    ...
  Stacktrace:
   [1] one(::Type{Any}) at ./missing.jl:105
   [2] reduce_empty(::typeof(*), ::Type{T} where T) at ./reduce.jl:308
   [3] reduce_empty(::typeof(Base.mul_prod), ::Type{T} where T) at ./reduce.jl:316
   [4] mapreduce_empty(::typeof(identity), ::Function, ::Type{T} where T) at ./reduce.jl:335
   [5] _mapreduce(::typeof(identity), ::typeof(Base.mul_prod), ::IndexLinear, ::Array{Any,1}) at ./reduce.jl:392
   [6] _mapreduce_dim at ./reducedim.jl:312 [inlined]
   [7] #mapreduce#580 at ./reducedim.jl:307 [inlined]
   [8] mapreduce at ./reducedim.jl:307 [inlined]
   [9] _prod at ./reducedim.jl:657 [inlined]
   [10] _prod(::Array{Any,1}, ::Colon) at ./reducedim.jl:656
   [11] #prod#585 at ./reducedim.jl:652 [inlined]
   [12] prod at ./reducedim.jl:652 [inlined]
   [13] parityinterior(::SubManifold{⟨+⟩,1,0x0000000000000001}, ::UInt64, ::UInt64) at /Users/eschnett/.julia/dev/Grassmann/src/parity.jl:77
   [14] interior(::UInt64, ::UInt64, ::SubManifold{⟨+⟩,1,0x0000000000000001}) at /Users/eschnett/.julia/dev/Grassmann/src/parity.jl:168
   [15] skewaddblade!_pre(::SubManifold{⟨+⟩,1,0x0000000000000001}, ::SizedArray{Tuple{1},Any,1,1}, ::UInt64, ::UInt64, ::Expr) at /Users/eschnett/.julia/dev/Grassmann/src/algebra.jl:258
   [16] #s335#570 at /Users/eschnett/.julia/dev/Grassmann/src/products.jl:373 [inlined]
   [17] #s335#570(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) at ./none:0
   [18] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any,N} where N) at ./boot.jl:526
   [19] dot(::Chain{v₁,1,Chain{v₁,1,Rational{Int128},1},1}, ::Chain{v₁,1,Rational{Int128},1}) at /Users/eschnett/.julia/packages/AbstractTensors/6BQWx/src/AbstractTensors.jl:132
   [20] top-level scope at /Users/eschnett/src/jl/RayTraceGR/test/runtests.jl:64
   [21] top-level scope at /Users/eschnett/src/julia-1.4/usr/share/julia/stdlib/v1.4/Test/src/Test.jl:1186
   [22] include(::String) at ./client.jl:439
   [23] top-level scope at none:6
   [24] eval(::Module, ::Any) at ./boot.jl:331
   [25] exec_options(::Base.JLOptions) at ./client.jl:264
   [26] _start() at ./client.jl:484

Tracking things down, the problem comes from an empty list of indices used in parity.jl line 77:

    ind = indices(B,N); g = prod(V[ind])

V[ind] is an empty array of type Vector{Any}, so prod cannot determine the neutral element. (I am using 1-dimensional manifolds, which is probably a bit of a special case.)

I believe the solution is to add explicit Int declarations to the array comprehensions in lines 198, 199 of DirectSum.jl in DirectSum:

@inline getindex(vs::SubManifold,i::Vector) = [getindex(vs,j) for j ∈ i]
@inline getindex(vs::SubManifold,i::UnitRange{Int}) = [getindex(vs,j) for j ∈ i]

This change resolves the problem for me.

[chakravala]

The type should not be explicitly declared, since it isn't always going to be Int. If you check with @code_warntype the method is already type stable, and sometimes it will not return an Int. For example, one might define the metric with Float64 instead of Int.

I will have to investigate this issue more, could you provide a minimum reproducible example>

[eschnett]

Found it... Was my mistake after all. The code below reproduces the problem:

using ComputedFieldTypes
using StaticArrays
using Grassmann

# @basis 1
V = v1
Vec(X) = Chain{V,1,X}
Mat(X) = Chain{V,1, fulltype(Vec(X))}
function Base.rand(::Type{<:Chain{V,1,T}}) where {V, T}
    D = ndims(V)
    Chain{V,1}(SVector{D,T}((rand(T) for i in 1:ndims(V))...))::Chain{V,1,<:T}
end
A = rand(Mat(Float64))
x = rand(Vec(Float64))

y = A ⋅ x

Of course, I didn't intend to define V in this way. If you remove this definition and use the @basis instead, the problem goes away.