Code Monkey home page Code Monkey logo

Comments (8)

MikaelSlevinsky avatar MikaelSlevinsky commented on June 20, 2024

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>

from fasttransforms.jl.

MikaelSlevinsky avatar MikaelSlevinsky commented on June 20, 2024

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>

from fasttransforms.jl.

Luapulu avatar Luapulu commented on June 20, 2024

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.

from fasttransforms.jl.

dlfivefifty avatar dlfivefifty commented on June 20, 2024

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

from fasttransforms.jl.

Luapulu avatar Luapulu commented on June 20, 2024

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.

from fasttransforms.jl.

dlfivefifty avatar dlfivefifty commented on June 20, 2024

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

from fasttransforms.jl.

Luapulu avatar Luapulu commented on June 20, 2024

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?

from fasttransforms.jl.

Luapulu avatar Luapulu commented on June 20, 2024

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.

from fasttransforms.jl.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.