juliaapproximation / fasttransforms.jl Goto Github PK

jjt([1.],1.5,-0.5,1.5,-0.5)

throw

LoadError: LoadError: UndefRefError: access to undefined reference

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

LowRankApprox.jl has a very large compile time and memory usage, due to the many mul! overrides involving StridedArrays.

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.

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.

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.

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.

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(n²).

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)

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

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.

@testset is preferable to println

There is an implementation here:

https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/src/SphericalHarmonics/Butterfly.jl#L121

but I feel like this could be improved.

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.

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).

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.

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?

Github: https://github.com/mweastwood/LibHealpix.jl

Docs: http://mweastwood.info/LibHealpix.jl/stable/

cjt, icjt, and jjt fail for coefficient vectors of length <= 2.

Complex and matrix-valued methods are missing.

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

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

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.

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.

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

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.

• Chebyshev T and U
• Jacobi (commit 65ffcdb)
• Spherical harmonics
• Triangle harmonics
• nfft
• Padua

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?

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?

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.

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 **

So that we can start using it in ApproxFun

femtocleaner https://github.com/JuliaComputing/FemtoCleaner.jl updates deprecated syntax automatically via a pull request.

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

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

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

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.

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.