Have `paduapts` return the data as the transpose of the current format?

[dlfivefifty]

Right now paduapts(n) returns an N x 2 matrix. I'm proposing that it instead return a 2 x N matrix. The reason is that this allows for a no-op conversion to a Vector{Vec{2,Float64}}: the current code works

pts=paduapts(10)
ptr=reinterpret(Ptr{Vec{2,Float64}},pointer(pts'))  # reinterpret the memory of pts'
unsafe_wrap(Vector{Vec{2,Float64}},
            ptr,
            size(pts,1),true)  

With the proposed change, the pts' would become pts and no new memory would need to be allocated.

[MikaelSlevinsky]

This doesn't work every time:

Richards-MBP-2:~ Mikael$ julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.3-pre.0 (2017-12-18 07:11 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 93168a6826* (90 days old release-0.6)
|__/                   |  x86_64-apple-darwin17.2.0

julia> using FastTransforms, StaticArrays

julia> pts=paduapoints(10)
66×2 Array{Float64,2}:
  1.0        0.959493
  1.0        0.654861
  1.0        0.142315
  1.0       -0.415415
  1.0       -0.841254
  1.0       -1.0     
  0.951057   1.0     
  0.951057   0.841254
  0.951057   0.415415
  0.951057  -0.142315
  ⋮                  
 -0.951057  -0.142315
 -0.951057  -0.654861
 -0.951057  -0.959493
 -1.0        0.959493
 -1.0        0.654861
 -1.0        0.142315
 -1.0       -0.415415
 -1.0       -0.841254
 -1.0       -1.0     

julia> ptr=reinterpret(Ptr{SVector{2,Float64}},pointer(pts'))  # reinterpret the memory of pts'
Ptr{SVector{2,Float64}} @0x00000001299d0e00

julia> unsafe_wrap(Vector{SVector{2,Float64}},
                   ptr,
                   size(pts,1),true)
66-element Array{SVector{2,Float64},1}:
 [2.27902e-314, 2.29691e-314]
 [2.29752e-314, 2.27908e-314]
 [2.29752e-314, 2.27902e-314]
 [2.29752e-314, 2.45615e-314]
 [2.29752e-314, 2.27902e-314]
 [2.27902e-314, 2.45615e-314]
 [2.27902e-314, 2.45616e-314]
 [2.29752e-314, 2.27908e-314]
 [2.29752e-314, 2.45615e-314]
 [2.27902e-314, 2.2969e-314] 
 ⋮                           
 [0.0, 0.0]                  
 [0.0, 0.0]                  
 [2.29752e-314, 2.27902e-314]
 [2.29752e-314, 2.27902e-314]
 [2.29752e-314, 2.45598e-314]
 [2.27902e-314, 2.45616e-314]
 [2.27902e-314, 2.45619e-314]
 [-1.0, -0.841254]           
 [-1.0, -1.0]                

julia> pts=paduapoints(10)
66×2 Array{Float64,2}:
  1.0        0.959493
  1.0        0.654861
  1.0        0.142315
  1.0       -0.415415
  1.0       -0.841254
  1.0       -1.0     
  0.951057   1.0     
  0.951057   0.841254
  0.951057   0.415415
  0.951057  -0.142315
  ⋮                  
 -0.951057  -0.142315
 -0.951057  -0.654861
 -0.951057  -0.959493
 -1.0        0.959493
 -1.0        0.654861
 -1.0        0.142315
 -1.0       -0.415415
 -1.0       -0.841254
 -1.0       -1.0     

julia> ptr=reinterpret(Ptr{SVector{2,Float64}},pointer(pts'))  # reinterpret the memory of pts'
Ptr{SVector{2,Float64}} @0x0000000129311bb0

julia> unsafe_wrap(Vector{SVector{2,Float64}},
                   ptr,
                   size(pts,1),true)
66-element Array{SVector{2,Float64},1}:
 [1.0, 0.959493]       
 [1.0, 0.654861]       
 [1.0, 0.142315]       
 [1.0, -0.415415]      
 [1.0, -0.841254]      
 [1.0, -1.0]           
 [0.951057, 1.0]       
 [0.951057, 0.841254]  
 [0.951057, 0.415415]  
 [0.951057, -0.142315] 
 ⋮                     
 [-0.951057, -0.142315]
 [-0.951057, -0.654861]
 [-0.951057, -0.959493]
 [-1.0, 0.959493]      
 [-1.0, 0.654861]      
 [-1.0, 0.142315]      
 [-1.0, -0.415415]     
 [-1.0, -0.841254]     
 [-1.0, -1.0]          

julia> 

[MikaelSlevinsky]

Hmm, but it's unrelated to the question. If I transpose the points first, then it works on the first shot:

julia> using FastTransforms, StaticArrays

julia> pts=paduapoints(10)'
2×66 Array{Float64,2}:
 1.0       1.0       1.0        1.0        1.0        1.0  0.951057  0.951057  0.951057   0.951057   0.951057  …  -0.951057  -0.951057  -0.951057  -0.951057  -0.951057  -1.0       -1.0       -1.0       -1.0       -1.0       -1.0
 0.959493  0.654861  0.142315  -0.415415  -0.841254  -1.0  1.0       0.841254  0.415415  -0.142315  -0.654861      0.841254   0.415415  -0.142315  -0.654861  -0.959493   0.959493   0.654861   0.142315  -0.415415  -0.841254  -1.0

julia> ptr=reinterpret(Ptr{SVector{2,Float64}},pointer(pts))  # reinterpret the memory of pts'
Ptr{SVector{2,Float64}} @0x00000001261ca040

julia> unsafe_wrap(Vector{SVector{2,Float64}},
                   ptr,
                   size(pts,2),true)
66-element Array{SVector{2,Float64},1}:
 [1.0, 0.959493]       
 [1.0, 0.654861]       
 [1.0, 0.142315]       
 [1.0, -0.415415]      
 [1.0, -0.841254]      
 [1.0, -1.0]           
 [0.951057, 1.0]       
 [0.951057, 0.841254]  
 [0.951057, 0.415415]  
 [0.951057, -0.142315] 
 [0.951057, -0.654861] 
 [0.951057, -0.959493] 
 [0.809017, 0.959493]  
 [0.809017, 0.654861]  
 [0.809017, 0.142315]  
 [0.809017, -0.415415] 
 [0.809017, -0.841254] 
 [0.809017, -1.0]      
 [0.587785, 1.0]       
 [0.587785, 0.841254]  
 [0.587785, 0.415415]  
 [0.587785, -0.142315] 
 [0.587785, -0.654861] 
 [0.587785, -0.959493] 
 [0.309017, 0.959493]  
 [0.309017, 0.654861]  
 [0.309017, 0.142315]  
 [0.309017, -0.415415] 
 [0.309017, -0.841254] 
 ⋮                     
 [-0.309017, 0.654861] 
 [-0.309017, 0.142315] 
 [-0.309017, -0.415415]
 [-0.309017, -0.841254]
 [-0.309017, -1.0]     
 [-0.587785, 1.0]      
 [-0.587785, 0.841254] 
 [-0.587785, 0.415415] 
 [-0.587785, -0.142315]
 [-0.587785, -0.654861]
 [-0.587785, -0.959493]
 [-0.809017, 0.959493] 
 [-0.809017, 0.654861] 
 [-0.809017, 0.142315] 
 [-0.809017, -0.415415]
 [-0.809017, -0.841254]
 [-0.809017, -1.0]     
 [-0.951057, 1.0]      
 [-0.951057, 0.841254] 
 [-0.951057, 0.415415] 
 [-0.951057, -0.142315]
 [-0.951057, -0.654861]
 [-0.951057, -0.959493]
 [-1.0, 0.959493]      
 [-1.0, 0.654861]      
 [-1.0, 0.142315]      
 [-1.0, -0.415415]     
 [-1.0, -0.841254]     
 [-1.0, -1.0]          

julia>

[Luapulu]

For my purposes I need multi dimensional transforms. So, I want to be able to Approximate a function from 2D to ND and transform the whole thing. In that case columns representing each dimension make the most sense to me. But I’ve also gone back and forth on this three times. Ultimately the initial point storage really doesn’t matter much, since it has to be copied into the plan anyways.

[dlfivefifty]

How does that relate to Padua points? It sounds like you want the standard tensor transforms.

[Luapulu]

Yeah, it doesn't really relate. The relevant point for me (which I didn't mention) is that to me each column ought to represent one dimension and each row ought to represent a point, just because that's how I think about tables and matrices. Because Julia uses column major arrays one may be tempted to transpose the matrix but unless there's actually a significant performance gain, I would be hesitant.

For users for whom this is actually a performance bottleneck I would be more in favour of offering an iterator of points which users may collect in any way they like. The default would be to write them in a matrix as is currently done.

[dlfivefifty]

I don't think iterator is the right approach, better to be a lazy AbstractVector{SVector{2,T}}.

[Luapulu]

For Padua points this doesn't work so well because the way you make them is by iterating through. You could rewrite the generation such that it uses explicit indexing but it would be a bit ugly. An iterator would be easier code.

In any case, is the creation of the points actually a performance bottleneck for anyone?

[Luapulu]

Here's the code in FastTransforms

FastTransforms.jl/src/PaduaTransform.jl

Lines 207 to 228 in 8be10bf

	function paduapoints(::Type{T}, n::Integer) where T
	N=div((n+1)*(n+2),2)
	MM=Matrix{T}(undef,N,2)
	m=0
	delta=0
	NN=fld(n+2,2)
	@inbounds for k=n:-1:0
	if isodd(n)>0
	delta=mod(k,2)
	end
	@inbounds for j=NN+delta:-1:1
	m+=1
	MM[m,1]=sinpi(T(k)/n-T(0.5))
	if isodd(n-k)>0
	MM[m,2]=sinpi((2j-one(T))/(n+1)-T(0.5))
	else
	MM[m,2]=sinpi(T(2j-2)/(n+1)-T(0.5))
	end
	end
	end
	return MM
	end

and here's my code

https://github.com/JuliaApproximation/ChebyshevTransforms.jl/blob/e8a7eac7c55ef6ecce5e12fc297d0bef4eb905e1/src/ChebyshevTransforms.jl#L191-L208

So, yeah, you can probably invert that to index explicitly but it might be a bit ugly.