juliaapproximation / fasttransforms.jl Goto Github PK
View Code? Open in Web Editor NEW:rocket: Julia package for orthogonal polynomial transforms :snowboarder:
License: Other
:rocket: Julia package for orthogonal polynomial transforms :snowboarder:
License: Other
This would mirror the current ApproxFun
transforms interface, and could help avoiding issues with converting At_mul_B!(Y,A,X)
overrides to mul!(Y,A',X)
. (I could be wrong, but it's not clear to me that the authors were aware of how this change in Base could affect packages that use recursive linear algebra data structures.) @dlfivefifty, thoughts?
I was talking to @dlfivefifty over at ApproxFun.jl
where they have an implementation of the hermite transform using gausshermite quadrature. While this is clearly a reliable choice, is there any scope for implementing an O(N*logN)
(i.e. faster than O(N^2)
) method as described in
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2630232/#!po=88.5000
?
For band limited functions there is a factorization using recursion the relation, akin to FFT vs basic FT. I don't know it well enough to know how practical it is to implement. On the surface it appears an O(N*logN)
method may be possible for any orthogonal polynomials with two-term recursion...
EDIT: it seems that the exact two-term recursion factorisation is numerically unstable, and in practice the transform is implemented as an approximate Hermite-Newton-Cotes transform on a linear grid. Probably my question should be closed here.
A stable method in O(mn) operations is proposed in
B. Adcock and A. C. Hansen. Stable reconstructions in Hilbert spaces and the resolution of the Gibbs phenomenon. Appl. Comput. Harmon. Anal., 32:357–388, 2012
where m is the number of Fourier modes and n is the number of Legendre modes. m = O(n2).
With a fresh install of the released version on Julia 0.6.2:
julia> using FastTransforms
INFO: Precompiling module FastTransforms.
WARNING: --output requested, but no modules defined during run
WARNING: The call to compilecache failed to create a usable precompiled cache file for module FFTW. Got:
WARNING: Cache file "/Users/dpsanders/.julia/lib/v0.6/FFTW.ji" not found.
WARNING: eval from module Main to ToeplitzMatrices:
Expr(:call, Expr(:., :Base, :include_from_node1)::Any, "/Users/dpsanders/.julia/v0.6/FFTW/src/FFTW.jl")::Any
** incremental compilation may be broken for this module **
WARNING: --output requested, but no modules defined during run
WARNING: The call to compilecache failed to create a usable precompiled cache file for module FFTW. Got:
WARNING: Cache file "/Users/dpsanders/.julia/lib/v0.6/FFTW.ji" not found.
WARNING: eval from module Main to LowRankApprox:
Expr(:call, Expr(:., :Base, :include_from_node1)::Any, "/Users/dpsanders/.julia/v0.6/FFTW/src/FFTW.jl")::Any
** incremental compilation may be broken for this module **
It errors out due to some zero length arrays. Changing div
to fld
helps one bug, but 0 x 0 Toeplitz matrices are broken JuliaLinearAlgebra/ToeplitzMatrices.jl#22
Any interest in having a nonuniform fast Fourier transform (FFT) and inverse NFFT in FastTransforms.jl
? I have one implemented, but it is not based on asymptotic expansions.
(v0.7) pkg> generate ISOLATION
Generating project ISOLATION:
ISOLATION/Project.toml
ISOLATION/src/ISOLATION.jl
(v0.7) pkg> activate .
(isolation) pkg> dev https://github.com/MikaelSlevinsky/FastTransforms.jl
Cloning default registries into /Users/vincentcp/.julia/registries
Cloning registry General from "https://github.com/JuliaRegistries/General.git"
Updating registry at `~/.julia/registries/General`
Updating git-repo `https://github.com/JuliaRegistries/General.git`
Cloning git-repo `https://github.com/MikaelSlevinsky/FastTransforms.jl`
Updating git-repo `https://github.com/MikaelSlevinsky/FastTransforms.jl`
┌ Warning: package at /var/folders/65/l0ldtg0j4yx5m57c8ztcs5yh0000gn/T/tmpdodGuF will need to have a [Julia]Project.toml file in the future
└ @ Pkg.Types /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v0.7/Pkg/src/Types.jl:619
Resolving package versions...
Installed AbstractFFTs ───────── v0.3.2
Installed Conda ──────────────── v0.8.1
...
Installed DataStructures ─────── v0.9.0
Updating `Project.toml`
[no changes]
Updating `Manifest.toml`
[no changes]
Building Conda ───────────→ `~/.julia/packages/Conda/S0nV/deps/build.log`
Building FFTW ────────────→ `~/.julia/packages/FFTW/Kjb9/deps/build.log`
Building SpecialFunctions → `~/.julia/packages/SpecialFunctions/kvjJ/deps/build.log`
julia> using FastTransforms
[ Info: Precompiling module FastTransforms
...
ERROR: LoadError: LoadError: LoadError: UndefVarError: REDFT01 not defined
Stacktrace:
[1] top-level scope at none:0
[2] include at ./boot.jl:317 [inlined]
[3] include_relative(::Module, ::String) at ./loading.jl:1034
[4] include at ./sysimg.jl:29 [inlined]
[5] include(::String) at /Users/vincentcp/.julia/dev/FastTransforms/src/FastTransforms.jl:2
[6] top-level scope at none:0
[7] include at ./boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at ./loading.jl:1034
[9] include at ./sysimg.jl:29 [inlined]
[10] include(::String) at /Users/vincentcp/.julia/dev/FastTransforms/src/FastTransforms.jl:2
[11] top-level scope at none:0
[12] include at ./boot.jl:317 [inlined]
[13] include_relative(::Module, ::String) at ./loading.jl:1034
[14] include(::Module, ::String) at ./sysimg.jl:29
[15] top-level scope at none:0
[16] eval at ./boot.jl:319 [inlined]
[17] eval(::Expr) at ./client.jl:394
[18] top-level scope at ./none:3 [inlined]
[19] top-level scope at ./<missing>:0
in expression starting at /Users/vincentcp/.julia/dev/FastTransforms/src/SphericalHarmonics/synthesisanalysis.jl:55
in expression starting at /Users/vincentcp/.julia/dev/FastTransforms/src/SphericalHarmonics/SphericalHarmonics.jl:16
in expression starting at /Users/vincentcp/.julia/dev/FastTransforms/src/FastTransforms.jl:90
ERROR: Failed to precompile FastTransforms to /Users/vincentcp/.julia/compiled/v0.7/FastTransforms/5Lm8.ji.
The types for the fft
functions seem to be overly specific.
(It looks like they should work for any Real
.)
Is this deliberate or should I loosen them?
running the readme examples
Julia 1.1:
using FastTransforms
F = sphrandn(Float64, 1024, 1024);
@time G = sph2fourier(F; sketch = :none);
Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:05
61.230271 seconds (598.27 M allocations: 10.007 GiB, 10.16% gc time)
@time H = fourier2sph(G; sketch = :none);
Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:05
64.233187 seconds (593.01 M allocations: 9.955 GiB, 10.55% gc time)
Julia 1.0.3:
F = sphrandn(Float64, 1024, 1024);
@time G = sph2fourier(F; sketch = :none);
Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:04
45.251658 seconds (436.73 M allocations: 7.941 GiB, 7.81% gc time)
@time H = fourier2sph(G; sketch = :none);
Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:03
49.485671 seconds (627.87 M allocations: 10.231 GiB, 11.12% gc time)
When using the sphevaluate function for the first time in a session, the result is incorrect. However, every subsequent time the function is used, the correct result is returned.
Below is a snippet of code that reproduces this issue.
I am using Julia version 0.6.2.
julia> using FastTransforms
julia> Pkg.status("FastTransforms")
- FastTransforms 0.3.0
julia> sphevaluate(0,0,0,0)
0.0209479177387814 + 0.0im
julia> sphevaluate(0,0,0,0)
0.28209479177387814 + 0.0im
Some seismic applications make use of partial sums involving a signal and the DFT matrix with variable summation limits. One approach decomposes the partial DFT hierarchically into rectangular subblocks that are related to fast convolution. This might require a HierarchicalDFTMatrix
LowRankApprox.jl has a very large compile time and memory usage, due to the many mul!
overrides involving StridedArray
s.
This will surely be fixed in time, but in the meantime, how hard would it be to remove the dependence here?
Another option would be to rewrite LowRankApprox.jl using LazyArrays. This has proven very successful for BandedMatrices: as of yesterday compile is down from 30s / 7GB to 2s / 100MB. But I think it's best to minimise the investment into that design.
I think other people might find the transforms useful in C and Python code. I think it's easy to wrap the Julia code in a C library:
http://docs.julialang.org/en/release-0.4/manual/embedding/
Down the line (0.5?) it will be possible to compile Julia code to a library that doesn't depend on an Julia installation.
jjt([1.],1.5,-0.5,1.5,-0.5)
throw
LoadError: LoadError: UndefRefError: access to undefined reference
julia> using FastTransforms
julia> ccheb = [34.0; 48.0; 18.0]
3-element Array{Float64,1}:
34.0
48.0
18.0
julia> cleg = cheb2leg(ccheb)
3-element Array{Float64,1}:
28.0
48.0
24.0
julia> cleg_TH = FastTransforms.th_cheb2leg(ccheb)
3-element Array{Float64,1}:
28.0
47.99999999999999
23.999999999999993
julia> ccheb = BigFloat[34.0; 48.0; 18.0]
3-element Array{BigFloat,1}:
3.4e+01
4.8e+01
1.8e+01
julia> cleg_TH = FastTransforms.th_cheb2leg(ccheb)
3-element Array{BigFloat,1}:
2.8e+01
-4.800000000000000000000000000000000000000000000000000000000000000000000000000111e+01
-2.399999999999999999999999999999999999999999999999999999999999999999999999999972e+01
All of the plan
types in the package were created on a need-to-plan basis. Now that things are coming together, it might be sensible to create a more effective type hierarchy.
FastTransforms.icjt(Float64[],0.0,0.0)
FMM acceleration of the Chebyshev--Legendre transform leads to a quasi-linear matrix-vector product. If an array of coefficients requires transformation, then it gets harder to beat level 3 BLAS because the dimensions of interest are square-rooted, due to memory limitations.
Currently, the data structures allow for an allocation-free implementation of matrix-vector products, since, for example, the LowRankMatrix
type includes a temporary vector variable re-used in intermediate steps. Perhaps a PlanMatrixMatrixMultiply
could allocate the necessary temporary arrays for the subblocks?
The butterfly algorithm is more complicated because the ranks of each subblock are determined adaptively.
@testset
is preferable to println
JuliaMath/AbstractFFTs.jl#28
(not sure which would be the better place to post this issue)
julia> VERSION # v"1.3.0-DEV.263"
julia> using DoubleFloats, FastTransforms, FFTW
julia> let df64s = rand(Double64, 2^10)
ComplexDF32( sum(df64s .- ifft(fft(df64s))) )
end
3.9988234e-30 + 0.0im
julia> VERSION # v"1.1.1"
julia> using DoubleFloats, FastTransforms, FFTW
julia> let df64s = rand(Double64, 2^10)
ComplexDF32( sum(df64s .- ifft(fft(df64s))) )
end
ERROR: type DoubleFloat{Float64} not supported
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] _fftfloat(::Type{DoubleFloat{Float64}}) at C:\Users\jas\.julia\packages\AbstractFFTs\PUqOK\src\definitions.jl:22
[3] _fftfloat(::DoubleFloat{Float64}) at C:\Users\jas\.julia\packages\AbstractFFTs\PUqOK\src\definitions.jl:23
[4] fftfloat(::DoubleFloat{Float64}) at C:\Users\jas\.julia\packages\AbstractFFTs\PUqOK\src\definitions.jl:18
[5] complexfloat(::Array{DoubleFloat{Float64},1}) at C:\Users\jas\.julia\packages\AbstractFFTs\PUqOK\src\definitions.jl:31
[6] fft(::Array{DoubleFloat{Float64},1}, ::UnitRange{Int64}) at C:\Users\jas\.julia\packages\AbstractFFTs\PUqOK\src\definitions.jl:198 (repeats 2 times)
[7] top-level scope at REPL[3]:2
The docs suggest that sph2fourier
and fourier2sph
are for complex-valued bases, but it looks like it's actual real-valued bases. This was our quick hack to convert to complex-valued SH:
function function2sph(f, n)
n = 20
θ = (0.5:n-0.5)/n * π
φ = (0:2n-2)*2/(2n-1) * π
F = [f(θ,φ) for θ in θ, φ in φ]
V = zero(F)
A_mul_B!(V, FastTransforms.plan_analysis(F), F)
M = size(V,2)
P = eye(Complex128,M)
for k = 2:4:M
P[k:k+1,k:k+1] = [im im;
-1 1]/sqrt(2)
end
# fix weird sign
for k = 4:4:M
P[k:k+1,k:k+1] = [-im im;
1 1]/sqrt(2)
end
(V |> fourier2sph) * P
end
I'm always confused which points to use for which transform, especially with the SphericalHarmonic/TriangleHarmonics case. I think an examples folder would help here.
If I figure out TriangleHarmonics, I'll add the example.
So that we can start using it in ApproxFun
According to R. Piessens and M. Branders, On the computation of Fourier transforms of singular functions, J. Comp. Appl. Math., 43:159-169, 1992, the modified Chebyshev moments of the complex exponential (multiplied by algebraic endpoint singularities) satisfy a many-term recurrence relation that is neither stable in the forward nor backward directions. How to proceed?
Combination of moment-based quadrature with fast orthogonal polynomial transforms to change bases is now irresistible: for example, the modified Legendre moments of the complex exponential are proportional to spherical Bessel functions that satisfy a three-term recurrence relation with well-known stability properties. By changing bases, we would simplify the recurrence and secure the stability by design.
julia> FastTransforms.cjt([1f0,2f0],0.0,0.0)
ERROR: MethodError: no method matching ForwardChebyshevUltrasphericalPlan(::Array{Float32,1}, ::Float64, ::Int64)
Closest candidates are:
ForwardChebyshevUltrasphericalPlan(::AbstractArray{T,1}, ::T, ::Int64) where T at /Users/solver/.julia/v0.6/FastTransforms/src/ChebyshevUltrasphericalPlan.jl:116
Stacktrace:
[1] (::FastTransforms.#kw##plan_cjt)(::Array{Any,1}, ::FastTransforms.#plan_cjt, ::Array{Float32,1}, ::Float64) at ./<missing>:0
[2] #plan_cjt#55 at /Users/solver/.julia/v0.6/FastTransforms/src/cjt.jl:164 [inlined]
[3] plan_cjt(::Array{Float32,1}, ::Float64, ::Float64) at /Users/solver/.julia/v0.6/FastTransforms/src/cjt.jl:164
[4] cjt(::Array{Float32,1}, ::Float64, ::Float64) at /Users/solver/.julia/v0.6/FastTransforms/src/cjt.jl:123
This issue is a placeholder for spherical harmonic bugs/enhancements that will be part of a patch:
the way the butterfly algorithm is currently implemented, it only works for dimensions that are powers of two, though nothing stops this from working for general dimensions. This afflicts FastSphericalHarmonicPlan
and ThinSphericalHarmonicPlan
. For the time being, it's recommended to just pad with zeros to the nextpow2
.
the commands sph2fourier
and fourier2sph
should work on different band limits in l
and m
.
the command plan_sph2fourier
is not type-stable as it very crudely dispatches to a different SphericalHarmonicPlan
based on the bandlimit. This is negligible compared to pre-computation costs for large degrees (and is somewhat likened to Base's factorize
).
the convenience constructors sphrand
, sphrandn
, etc... take m
and n
and create an array with dimensions m
and 2n-1
which is arguably misleading.
it would be great to have an allocation-free FFTW
override for converting the bivariate Fourier array to function values on the sphere at tensor product grids equispaced in angle. (Added starting with 67a6b56)
remove segmentation fault when using an even number of columns, thanks to @meggart for finding this. (Added starting with c0e4a54)
add methods for complex data,
and, NUFFT
variants for nonuniform point.
cjt
, icjt
, and jjt
fail for coefficient vectors of length <= 2.
Complex and matrix-valued methods are missing.
And if so, what kind of fast transform is applicable? I've heard they're related to Jacobi polynomials with negative parameters, but they might transform nicely to Taylor polynomials.
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.
The FFT provided in FastTransforms is specific to BigFloat
. It would be nice to support other abstract float types. Any reason why AbstractFloats
is not defined like this:
const AbstractFloats = Union{AbstractFloat,Complex{<:AbstractFloat}}
?
I ran a test using the DoubleFloats.jl
package. There are a few places where BigFloat
is hardcoded in genericfft.jl
, but apart from that the only issues I encountered were associated either with DoubleFloats
or Julia Base (see JuliaMath/DoubleFloats.jl#65).
julia> cjt([1.,2],.5,.5)
2-element Array{Float64,1}:
NaN
3.0
The Toeplitz-dot-Hankel approach has now been added and it's blazingly fast. But with two types of low-precomputation fast transforms, we can now run comparisons to extremely high order.
This gist uses normally distributed coefficients (seeded for reproducibility) and compares absolute error in the coefficients of the compositions of forward-inverse and inverse-forward Chebyshev--Legendre transforms.
The culprit appears to be the recently added cheb2leg
:
julia> using FastTransforms
julia> srand(0);
julia> c = randn(100_000);
julia> norm(cjt(c,0.,0.)-leg2cheb(c),Inf)
8.992806499463768e-15
julia> norm(icjt(c,0.,0.)-cheb2leg(c),Inf)
1.2517668892542133e-6
Is there a transform from a grid to the basis expected in fourier2ph
?
I can use 2D Fourier and drop coefficients, but that seems wasteful.
There is an implementation here:
but I feel like this could be improved.
Hi Mikael,
I am trying to understand how difficult it would be to use FastTransforms.jl to compute spherical harmonic transforms of healpix maps. Healpix maps have pixels arranged along rings of constant latitude, which the C++ Healpix library uses to give me "fast" spherical harmonic transforms. I would really like to compare your spherical harmonic transform against the existing one, because I anticipate yours is probably much faster.
Unfortunately I'm having a hard time understanding what format your routines need. Are there any plans to expand and document the synthesis and analysis code? Can it handle the case where there are a different number of pixels on each ring?
femtocleaner https://github.com/JuliaComputing/FemtoCleaner.jl updates deprecated syntax automatically via a pull request.
On my 32 bit linux machine, the precompilation failed under Julia v1.0.1 as follows. On the other hand, I didn't have any problems on 64 bit linux machine.
julia> using FastTransforms
[ Info: Precompiling FastTransforms [057dd010-8810-581a-b7be-e3fc3b93f78c]
ERROR: Failed to precompile FastTransforms [057dd010-8810-581a-b7be-e3fc3b93f78c] to /home/xxx/.julia/compiled/v1.0/FastTransforms/5Lm8s.ji.
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] macro expansion at ./logging.jl:313 [inlined]
[3] compilecache(::Base.PkgId, ::String) at ./loading.jl:1187
[4] _require(::Base.PkgId) at ./logging.jl:311
[5] require(::Base.PkgId) at ./loading.jl:855
[6] require(::Module, ::Symbol) at ./logging.jl:311
Currently I don't see a way to compute the nonuniform fourier transform at arbitrary points. nufft2
is the closest variant, but it requires a regular grid of x, y
coordinates in the transform space. A more flexible one should take any number of arbitrary (x, y)
pairs and return the same number of fourier coefficients. Or did I overlook something?
Amazing package @MikaelSlevinsky!
I'm wondering if it is possible to add a plan_fourier2sph
function. Also it would be helpful if you had an example of how to use current plan_sph2fourier
.
F = zeros(10,10)
F[1,1] = 1
F̃ = cheb2tri(F, -0.0, -0.5, -0.5)
F̃[1,1] ≈ sqrt(2π) # My oh my...why?
This looks like a really great package. You have a great README and tests - I was wondering if you upload coverage statistics and if so, would you be willing to put a badge on your README?
Cheers,
Katharine
julia> s = rand(BigFloat, 10);
julia> conv(s,s)
ERROR: BoundsError: attempt to access 108-element Array{BigFloat,1} at index [109]
julia> size(s)
(54,)
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.