Type inference fails for `GaussLegendre` integrals

[mikeingold]

I just used Cthulhu.descend to quickly peek through the package and noticed that integral(f, ::Segment, ::GaussLegendre) inferences an output type ::Any. The _integral(..., ::GaussLegendre) function looks almost entirely blue (good), except for this line

MeshIntegrals.jl/src/integral.jl

Line 93 in c7c0a47

return (1 // (2^N)) .* sum(integrand, zip(weight_grid, node_grid))

where I see an inference zip(weight_grid, node_grid)::Any. However, I can't descend any further into this line because I get a This is a runtime call notice for the line itself and nothing beneath zip.

My experience deep-diving type inference is very limited, so I'm not sure if this is much of an issue, or how to go about addressing it beyond what I've already done.

Details

Cthulhu.jl output:

68 function (_integral(
69         f::Core.Const(Main.f),
70         geometry::Segment{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}},
71         rule::GaussLegendre::GaussLegendre;
72         FP::Type{Float64}::Type{T} = Float64,
73         diff_method::FiniteDifference{Float64}::DM = _default_diff_method(geometry)
74 ) where {DM <: DifferentiationMethod, T <: AbstractFloat})::Any
75     N::Core.Const(1) = Meshes.paramdim(geometry)
77
77     # Get Gauss-Legendre nodes and weights of type FP for a region [-1,1]ᴺ
78     xs = Iterators.map(FP::Type{Float64}, rule::GaussLegendre.nodes::Vector{Float64})
79     ws = Iterators.map(FP::Type{Float64}, rule::GaussLegendre.weights::Vector{Float64})
80     weight_grid = Iterators.product(ntuple(Returns(ws), N::Core.Const(1))...)
81     node_grid = Iterators.product(ntuple(Returns(xs), N::Core.Const(1))...)
83
83     # Domain transformation: x [-1,1] ↦ t [0,1]
84     t(x) = (1 // 2) * x + (1 // 2)
86
86     function integrand((weights, nodes))
87         # ts = t.(nodes), but non-allocating
88         ts = ntuple(i -> t(nodes[i]), length(nodes))
89         # Integrand function
90         prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)
91     end
93
93     return (((1 // ((2^N::Core.Const(1))::Core.Const(2)))::Core.Const(1//2)) .*::Core.Const(*) sum(integrand::MeshIntegrals.var"#integrand#28"{FiniteDifference{Float64}, typeof(f), Segment{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, MeshIntegrals.var"#t#27"}, zip(weight_grid, node_grid))::Any)::Any