Extension of example from FAQ "interpolate a function at a specified grid" to 3D

[CalebBell]

Hello Dr. Olver, I've been seeking to use ApproxFun to approximate a variety of functions and had a few questions. My functions are expensive to compute, so I need to evaluate them on a fixed number of points evaluated once only and then run the fitting. I've been seeking to use the fixed number of points variant as described in the FAQ, which gives a working example of a 2D approximation based computing a Vandermonde matrix. I need to create approximations in 3D however. The extension to 3D would seem to be straightforward, but I'm running into errors. Is this type of approximation within the capabilities of ApproxFun at present?

S = Chebyshev(0..1)^3;
n = 1000; m = 50;
srand(0); x = rand(n); y = rand(n); z = rand(n);
v = exp.(x .* cos.(y) .* sin.(z));  # values at the non-default grid

V = Array{Float64}(n,m); # Create a Vandermonde matrix by evaluating the basis at the grid

for k = 1:m
    V[:,k] = Fun(S, [zeros(k-1);1]).(x, y, z)
end

f = Fun(S, V\v);

The error message is as follows:

MethodError: no method matching cumsum(::ApproxFun.CumSumIterator{ApproxFun.UnitCount{Int64}})
Closest candidates are:
  cumsum(::ApproxFun.Repeated{Bool}) at /home/caleb/.julia/v0.6/ApproxFun/src/LinearAlgebra/helper.jl:1163
  cumsum(::AbstractArray{T,N} where N) where T at multidimensional.jl:611
  cumsum(::AbstractArray{T,N} where N, ::Integer) where T at multidimensional.jl:611
  ...

Stacktrace:
 [1] block(::ApproxFun.Tensorizer{Tuple{ApproxFun.Repeated{Bool},ApproxFun.Repeated{Bool},ApproxFun.Repeated{Bool}}}, ::Int64) at /home/caleb/.julia/v0.6/ApproxFun/src/Multivariate/TensorSpace.jl:106
 [2] totensor(::ApproxFun.Tensorizer{Tuple{ApproxFun.Repeated{Bool},ApproxFun.Repeated{Bool},ApproxFun.Repeated{Bool}}}, ::Array{Float64,1}) at /home/caleb/.julia/v0.6/ApproxFun/src/Multivariate/TensorSpace.jl:566
 [3] evaluate(::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64}},ApproxFun.ProductDomain{Tuple{ApproxFun.Segment{Float64},ApproxFun.Segment{Float64},ApproxFun.Segment{Float64}},SVector{3,Float64}},Float64},Float64,Array{Float64,1}}, ::Float64, ::Float64, ::Float64) at /home/caleb/.julia/v0.6/ApproxFun/src/Fun/Fun.jl:189
 [4] (::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64}},ApproxFun.ProductDomain{Tuple{ApproxFun.Segment{Float64},ApproxFun.Segment{Float64},ApproxFun.Segment{Float64}},SVector{3,Float64}},Float64},Float64,Array{Float64,1}})(::Float64, ::Vararg{Float64,N} where N) at /home/caleb/.julia/v0.6/ApproxFun/src/Fun/Fun.jl:192
 [5] broadcast_t(::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64}},ApproxFun.ProductDomain{Tuple{ApproxFun.Segment{Float64},ApproxFun.Segment{Float64},ApproxFun.Segment{Float64}},SVector{3,Float64}},Float64},Float64,Array{Float64,1}}, ::Type{Any}, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N} where N) at ./broadcast.jl:256
 [6] broadcast_c(::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64}},ApproxFun.ProductDomain{Tuple{ApproxFun.Segment{Float64},ApproxFun.Segment{Float64},ApproxFun.Segment{Float64}},SVector{3,Float64}},Float64},Float64,Array{Float64,1}}, ::Type{Array}, ::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N} where N) at ./broadcast.jl:319
 [7] broadcast(::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64}},ApproxFun.ProductDomain{Tuple{ApproxFun.Segment{Float64},ApproxFun.Segment{Float64},ApproxFun.Segment{Float64}},SVector{3,Float64}},Float64},Float64,Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N} where N) at ./broadcast.jl:434
 [8] macro expansion at ./In[90]:9 [inlined]
 [9] anonymous at ./<missing>:?
 [10] include_string(::String, ::String) at ./loading.jl:515

I'm running Julia 0.6 and ApproxFun 0.7.0. Any help would be greatly appreciated. I've been incredibly impressed with the 1D capabilities of ApproxFun and used them on several occasions. Rational approximation is incredibly powerful; there are many problems in Chemical Engineering I hope to apply ApproxFun to! Thank you for contributing this excellent project!

Sincerely, Caleb

[dlfivefifty]

ApproxFun doesn't really support 3D yet, but I'll try to get your example working.

[CalebBell]

Hello Dr. Olver, Thank you for the kind reply. I would be very appreciative of the feature.
Please don't reprioritize this on my behalf though! ApproxFun is an incredibly large undertaking, and I'm sure you know know best where to lead it, as your time and interest permits. Sharing your excellent 1D and 2D support is incredibly kind already, and it is a testament to your APIs clear design that I thought ApproxFun might support 3D approximation too. No worries at all. Thanks again, Caleb

[dlfivefifty]

I had a look and it's going go be harder to fix than I thought: though the user facing code is indeed written in a way to be extendable to 3D (and higher), most of the implementations (like evaluate) are hard-coded for 2D.

There's some planned revamps on how coefficients are stored that will make this task easier, so I think it'll have to wait for now.

Glad the 1D and 2D code has been useful, though!