FFTs with posits #64

milankl · 2022-05-07T13:52:57Z

Base.maxintfloat(::Type{Posit}) has to be defined.

julia> using SoftPosit, FastTransforms

julia> x = Posit16.(rand(8))
8-element Vector{Posit16}:
 Posit16(0.6254883)
 Posit16(0.7780762)
 Posit16(0.8145752)
 Posit16(0.50598145)
 Posit16(0.74768066)
 Posit16(0.29797363)
 Posit16(0.89624023)
 Posit16(0.6437988)

julia> fft(x)
ERROR: MethodError: no method matching maxintfloat(::Type{Posit16})

The text was updated successfully, but these errors were encountered:

milankl · 2022-05-10T13:24:44Z

I believe maxintfloats are

Base.maxintfloat(::Type{Posit8}) = 8              # = 2^3
Base.maxintfloat(::Type{Posit16}) = 512           # = 2^9
Base.maxintfloat(::Type{Posit16_2}) = 1024        # = 2^10
Base.maxintfloat(::Type{Posit32}) = 8388608       # = 2^23

milankl · 2022-06-10T17:08:12Z

Base.maxintfloat(::AbstractPosit) added with #68

milankl · 2022-06-10T17:35:36Z

@tkgunaratne let me know if there's anything else needed that we can include in v0.5 for FFTs with posits

tkgunaratne · 2022-06-10T19:27:08Z

success!!

Thanks a lot!! :-)

giordano · 2022-06-10T19:47:23Z

Looks like they're converted to Complex{Float32}?

tkgunaratne · 2022-06-10T19:50:46Z

Oops didn't notice that :-(. Got excited as the FFT function worked.

milankl · 2022-06-11T14:09:04Z

That must be an explicit conversion as I deliberately do not define promotion between floats and posits

julia> 1f0*Posit16(1)
ERROR: promotion of types Float32 and Posit16 failed to change any arguments

milankl · 2022-06-11T14:21:35Z

sinpi is missing

Base.sinpi(x::T) where {T<:AbstractPosit} = convert(T,sinpi(float(x)))

then fft still seems to convert to Float32, but calling generic_fft_pow2 does the job

julia> using FastTransforms

julia> x = complex(Posit16.(randn(8)));

julia> FastTransforms.generic_fft_pow2(x)
8-element Vector{Complex{Posit16}}:
  Posit16(3.8544922) + Posit16(0.0)im
  Posit16(2.0048828) - Posit16(2.7148438)im
 Posit16(-3.2929688) + Posit16(0.78149414)im
  Posit16(1.4091797) - Posit16(1.3051758)im
  Posit16(3.6767578) + Posit16(0.0)im
  Posit16(1.4082031) + Posit16(1.3061523)im
 Posit16(-3.2929688) - Posit16(0.78149414)im
  Posit16(2.0039062) + Posit16(2.7148438)im

giordano · 2022-06-11T14:28:57Z

As a sanity check I usually do the round-trip fft->ifft:

julia> using FastTransforms, SoftPosit

julia> Base.sinpi(x::T) where {T<:AbstractPosit} = convert(T,sinpi(float(x)))

julia> x = complex(Posit16.(randn(8)))
8-element Vector{Complex{Posit16}}:
 Posit16(-0.9309082) + Posit16(0.0)im
 Posit16(0.44580078) + Posit16(0.0)im
 Posit16(-2.0751953) + Posit16(0.0)im
 Posit16(-1.7192383) + Posit16(0.0)im
 Posit16(-2.0234375) + Posit16(0.0)im
 Posit16(-0.2477417) + Posit16(0.0)im
 Posit16(-1.5751953) + Posit16(0.0)im
 Posit16(-1.5771484) + Posit16(0.0)im

julia> FastTransforms.generic_ifft_pow2(FastTransforms.generic_fft_pow2(x))
8-element Vector{Complex{Posit16}}:
 Posit16(-0.93115234) + Posit16(0.0)im
  Posit16(0.44604492) - Posit16(0.0002746582)im
  Posit16(-2.0742188) + Posit16(0.0)im
  Posit16(-1.7197266) + Posit16(0.00022888184)im
  Posit16(-2.0234375) + Posit16(0.0)im
 Posit16(-0.24731445) - Posit16(0.00015258789)im
  Posit16(-1.5751953) + Posit16(0.0)im
  Posit16(-1.5771484) + Posit16(0.00019836426)im

Look close enough. Note that isapprox doesn't work on Posit16 because it doesn't define eps.

milankl · 2022-06-13T12:54:22Z

Forgot to add Base. ...

SoftPosit.jl/src/arithmetics.jl

Lines 45 to 49 in 7f32549

    
           eps(::Type{Posit8}) = reinterpret(Posit8,0x28) 
        
           eps(::Type{Posit16}) = reinterpret(Posit16,0x0a00) 
        
           eps(::Type{Posit16_1}) = reinterpret(Posit16_1,0x0100) 
        
           eps(::Type{Posit32}) = reinterpret(Posit32,0x00a0_0000) 
        
           eps(x::AbstractPosit) = max(x-prevfloat(x),nextfloat(x)-x)

milankl · 2022-06-13T14:02:35Z

We still can't do isapprox(::Vector{Posit16},Vector{Posit16}) because

julia> using .SoftPosit
julia> using FastTransforms
julia> x = complex(Posit16.(randn(8)))
8-element Vector{Complex{Posit16}}:
   Posit16(2.196289) + Posit16(0.0)im
 Posit16(0.31066895) + Posit16(0.0)im
  Posit16(0.6977539) + Posit16(0.0)im
  Posit16(1.1914062) + Posit16(0.0)im
  Posit16(0.2701416) + Posit16(0.0)im
 Posit16(-0.9902344) + Posit16(0.0)im
 Posit16(0.47131348) + Posit16(0.0)im
  Posit16(0.2368164) + Posit16(0.0)im

julia> x2 = FastTransforms.generic_ifft_pow2(FastTransforms.generic_fft_pow2(x))
8-element Vector{Complex{Posit16}}:
   Posit16(2.1972656) + Posit16(0.0)im
  Posit16(0.31079102) + Posit16(6.67572e-6)im
   Posit16(0.6977539) + Posit16(0.0)im
   Posit16(1.1914062) - Posit16(0.00025081635)im
  Posit16(0.27026367) + Posit16(0.0)im
 Posit16(-0.98999023) + Posit16(0.00025081635)im
   Posit16(0.4711914) + Posit16(0.0)im
  Posit16(0.23693848) - Posit16(6.67572e-6)im

julia> x ≈ x2
ERROR: promotion of types Float64 and Posit16 failed to change any arguments
Stacktrace:
  [1] error(::String, ::String, ::String)
    @ Base ./error.jl:42
  [2] sametype_error(input::Tuple{Float64, Posit16})
    @ Base ./promotion.jl:374
  [3] not_sametype(x::Tuple{Float64, Posit16}, y::Tuple{Float64, Posit16})
    @ Base ./promotion.jl:368
  [4] promote
    @ ./promotion.jl:351 [inlined]
  [5] +(x::Float64, y::Posit16)
    @ Base ./promotion.jl:379
  [6] generic_norm2(x::Vector{Complex{Posit16}})
    @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:507

it seems that the LinearAlgebra.normInf(x::Posit16) returns a Float32 that gets promoted to Float64 here which causes a problem in an addition in LinearAlgebra.generic_norm2.

Otherwise isapprox works fine

julia> a = Posit16(rand())
julia> a ≈ a
true

julia> complex(a) ≈ complex(a)
true

tkgunaratne · 2022-06-13T14:37:10Z

Could this line "julia> Base.sinpi(x::T) where {T<:AbstractPosit} = convert(T,sinpi(float(x)))" be added to any file in SoftPosit/src?

giordano · 2022-06-13T14:37:47Z

It's already in #70 😉

tkgunaratne · 2022-06-22T18:39:26Z

I was just wondering whether any additional work has to be done to support for non-power of 2 fft/iffts with SoftPosit?

milankl · 2022-06-23T12:00:46Z

This https://github.com/JuliaDSP/FourierTransforms.jl looks like a generic not-just-power2 FFT implementation to me. It hasn't been touched in two years though. I think in general, we should probabaly try to push for a generic Julia FFT in a standalone package, as has been previously discussed in other places (JuliaApproximation/FastTransforms.jl#86). I'd be happy to help out, although as a first step I'd vote to make power2 FFTs fast before tackling non-power2s, but maybe that could all go hand in hand.

milankl · 2022-06-23T13:44:30Z

Having said that, we soon need a generic FFT for SpeedyWeather.jl which currently uses FFTW.rFFTWPlan{NF} and FFTW.rFFTWPlan{Complex{NF}} from FFTW.jl which only supports NF=Float32/64. But again the focus here would be rather to extend the FFTW interface to generic types than to support non-power2. If you are interested though, @tkgunaratne we could have a call and see what's actually needed to create such a package.

giordano · 2022-06-23T13:47:44Z

https://github.com/JuliaDSP/FourierTransforms.jl Edit: oh, it's the same one you mentioned above, but it has been touched just two weeks ago!

tkgunaratne · 2022-08-20T05:29:14Z

I'm afraid that I upgraded to Julia 1.8.0 and that seemed to have broken the fft operation with Posits :-(

ERROR: promotion of types Float64 and Posit16 failed to change any arguments Stacktrace: [1] error(::String, ::String, ::String) @ Base ./error.jl:44 [2] sametype_error(input::Tuple{Float64, Posit16}) @ Base ./promotion.jl:383 [3] not_sametype(x::Tuple{Float64, Posit16}, y::Tuple{Float64, Posit16}) @ Base ./promotion.jl:377 [4] promote @ ./promotion.jl:360 [inlined] [5] <=(x::Float64, y::Posit16) @ Base ./promotion.jl:429 [6] >=(x::Posit16, y::Float64) @ Base ./operators.jl:429 [7] sinpi(x::Posit16) @ Base.Math ./special/trig.jl:758 [8] generic_fft_pow2!(x::Vector{Posit16}) @ GenericFFT ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:106 [9] generic_fft_pow2(x::Vector{Complex{Posit16}}) @ GenericFFT ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:126 [10] generic_fft(x::Vector{Complex{Posit16}}) @ GenericFFT ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:50 [11] generic_fft @ ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:22 [inlined] [12] * @ ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:260 [inlined] [13] fft(x::Vector{Complex{Posit16}}) @ AbstractFFTs ~/.julia/packages/AbstractFFTs/Wg2Yf/src/definitions.jl:62 [14] top-level scope @ ~/Julia_Scripts/SoftPosit.jl/Examples/FFT_Example.jl:6

milankl · 2022-08-20T12:29:29Z

I cannot reproduce this issue

julia> using SoftPosit, GenericFFT

julia> fft(Posit16.(randn(8)))
8-element Vector{Complex{Posit16}}:
 Posit16(-0.24487305) + Posit16(0.0)im
   Posit16(4.1367188) - Posit16(0.7182617)im
   Posit16(0.8757324) - Posit16(1.0927734)im
  Posit16(-1.3613281) + Posit16(0.60791016)im
   Posit16(0.2409668) + Posit16(0.0)im
  Posit16(-1.3603516) - Posit16(0.6088867)im
   Posit16(0.8757324) + Posit16(1.0927734)im
    Posit16(4.138672) + Posit16(0.7192383)im

with Julia v1.8, SoftPosit v0.5, GenericFFT v0.1.1. Also it seems to have been triggered from >=(x::Posit16, y::Float64) in sinpi(x::Posit16), which in Julia v1.8 is fine

julia> sinpi(Posit16(1.5))
Posit16(-1.0)

Note in general that an error related to promotion is not due to SoftPosit.jl as we deliberately don't implement promotion between floats and posits to flag where unintentional type conversions happen (as users of this package are probably interested to do all arithmetic in posits and not floats).

tkgunaratne · 2022-10-11T08:25:11Z

I greatly appreciate your prompt response. It seems that I have been using Sofposit version 0.4 :-(. It seems to work fine now. I may take this opportunity to ask if there is any progress on non power of 2 point FFTs.

…

On Sat, Aug 20, 2022 at 5:29 AM Milan Klöwer ***@***.***> wrote: I cannot reproduce this issue julia> using SoftPosit, GenericFFT julia> fft(Posit16.(randn(8)))8-element Vector{Complex{Posit16}}: Posit16(-0.24487305) + Posit16(0.0)im Posit16(4.1367188) - Posit16(0.7182617)im Posit16(0.8757324) - Posit16(1.0927734)im Posit16(-1.3613281) + Posit16(0.60791016)im Posit16(0.2409668) + Posit16(0.0)im Posit16(-1.3603516) - Posit16(0.6088867)im Posit16(0.8757324) + Posit16(1.0927734)im Posit16(4.138672) + Posit16(0.7192383)im with Julia v1.8, SoftPosit v0.5, GenericFFT v0.1.1. Also it seems to have been triggered from >=(x::Posit16, y::Float64) in sinpi(x::Posit16), which in Julia v1.8 is fine julia> sinpi(Posit16(1.5))Posit16(-1.0) Note in general that an error related to promotion is not due to SoftPosit.jl as we deliberately don't implement promotion between floats and posits to flag where unintentional type conversions happen (as users of this package are probably interested to do all arithmetic in posits and not floats). — Reply to this email directly, view it on GitHub <#64 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AUFYAAJHOD3RRWLM5GPVGLDV2DFTJANCNFSM5VKOLPDA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

milankl · 2022-10-11T10:27:38Z

It seems that

julia> using SoftPosit, GenericFFT

julia> fft(Posit16.(randn(10)))
ERROR: MethodError: no method matching Posit16(::Irrational{:π})

is missing, which can quickly be added by conversion first to Float32/64

julia> Posit16(x::Irrational) = Posit16(Base.floattype(Posit16)(x))

But then there's also

julia> fft(Posit16.(randn(10)))
ERROR: promotion of types Float64 and Posit16 failed to change any arguments
Stacktrace:
  [1] error(::String, ::String, ::String)
    @ Base ./error.jl:44
  [2] sametype_error(input::Tuple{Float64, Posit16})
    @ Base ./promotion.jl:383
  [3] not_sametype(x::Tuple{Float64, Posit16}, y::Tuple{Float64, Posit16})
    @ Base ./promotion.jl:377
  [4] promote
    @ ./promotion.jl:360 [inlined]
  [5] *(x::Float64, y::Posit16)
    @ Base ./promotion.jl:389
  [6] lerpi
    @ ./range.jl:958 [inlined]
  [7] unsafe_getindex(r::LinRange{Posit16, Int64}, i::Int64)
    @ Base ./range.jl:951
  ...
 [19] getindex
    @ ./broadcast.jl:597 [inlined]
 [20] macro expansion
    @ ./broadcast.jl:961 [inlined]
 [21] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [22] copyto!
    @ ./broadcast.jl:960 [inlined]
 [23] copyto!
    @ ./broadcast.jl:913 [inlined]
 [24] copy
    @ ./broadcast.jl:885 [inlined]
 [25] materialize
    @ ./broadcast.jl:860 [inlined]
 [26] generic_fft(x::Vector{Complex{Posit16}})
    @ GenericFFT ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:52
 [27] generic_fft
    @ ~/.julia/packages/GenericFFT/Oa0Sx/src/fft.jl:22 [inlined]

Which has to be addressed in GenericFFT.jl. Non-power2 therein is implemented as a direct transform, so the complexity is worse than with FFTW, but it works (using Float16 which doesn't have the problems discussed above).

julia> GenericFFT.generic_fft(rand(Float16,8))
8-element Vector{ComplexF16}:
   Float16(5.895) + Float16(0.0)im
   Float16(-0.47) - Float16(0.3445)im
 Float16(-0.3438) - Float16(0.2461)im
  Float16(0.3042) + Float16(0.05566)im
 Float16(-0.7695) + Float16(0.0)im
  Float16(0.3044) - Float16(0.0569)im
 Float16(-0.3438) + Float16(0.2461)im
 Float16(-0.4688) + Float16(0.3457)im

julia> GenericFFT.generic_fft(rand(Float16,10))
10-element Vector{ComplexF16}:
    Float16(4.754) - Float16(0.001587)im
   Float16(0.6123) - Float16(0.1678)im
 Float16(-0.09766) + Float16(1.582)im
     Float16(0.28) - Float16(0.722)im
  Float16(-0.2092) - Float16(0.3606)im
     Float16(0.35) + Float16(0.00579)im
  Float16(-0.2206) + Float16(0.3496)im
   Float16(0.2737) + Float16(0.7065)im
 Float16(-0.05493) - Float16(1.603)im
    Float16(0.633) + Float16(0.1726)im

So if you are interested to get this going, we should consider updating GenericFFT.jl. Shouldn't be hard, the package probably just needs some care.

tkgunaratne · 2022-10-11T17:21:06Z

Appreciate it Milan!

However, when I tried julia> GenericFFT.generic_fft(rand(Float16,150)) I got the following

150-element Vector{ComplexF16}:  
 NaN16 + NaN16*im  
 NaN16 + NaN16*im  
 NaN16 + NaN16*im  
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im
       ⋮
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im
 NaN16 + NaN16*im

This opposed to

julia> GenericFFT.generic_fft(rand(Float16,256))
256-element Vector{ComplexF16}:
  Float16(128.2) + Float16(0.0)im
   Float16(3.02) + Float16(4.0)im
 Float16(-5.164) - Float16(2.139)im
 Float16(0.3098) - Float16(1.306)im
  Float16(1.049) + Float16(1.336)im
 Float16(-1.275) - Float16(4.324)im
  Float16(3.809) - Float16(1.603)im
                 ⋮
  Float16(3.797) + Float16(1.619)im
 Float16(-1.293) + Float16(4.324)im
  Float16(1.053) - Float16(1.318)im
 Float16(0.3127) + Float16(1.313)im
  Float16(-5.18) + Float16(2.129)im
  Float16(3.031) - Float16(3.996)im

Also, I noticed steep degradation of accuracy when compared with the Float64 implementation.

Also, I got following warning

┌ Warning: Using generic fft for FFTW number type.
 @ GenericFFT C:\Users\gunaratnet\.julia\packages\GenericFFT\Oa0Sx\src\fft.jl:48

milankl · 2022-10-13T12:42:35Z

Yes, don't forget that there's two things here. 1) Get a generic fft to return a Vector{Complex{T}} for any-length Vector{T} that's provided as argument and 2) address precision and range issues within that algorithm.

We have almost ticked off 1), but 2) requires more work. I'm not surprised that for large vectors the direct FT may hit overflows with Float16, which I believe is because it's a long sum of elements that may easily go beyond floatmax(Float16). However, that would be important to know, and clearly points towards the need for a recursive algorithm as the FFT is

tkgunaratne · 2022-10-13T16:57:32Z

I agree, I think there should be some modification done to the FFT algorithm to cater for Float16 Posit16 formats.

floatmax(Float16) = 65504 and max(abs(fft(rand(150)))) ~ 150

tkgunaratne · 2023-04-23T20:59:47Z

Getting back to the topic of non power of two FFTs;
MicroFloatingPoints.jl managed to facilitate non power of two FFTs for both GenericFFT.jl and FastTransforms.jl packages with few simple modifications to its code. I suspect it was fairly straight forward because under the hood it is all IEE 754 compatible Floats.

tkgunaratne · 2024-09-17T19:27:01Z

Hi @milankl, Just wondering whether it is too difficult to support non-power of 2 fft/iffts with Posits :-D.

milankl mentioned this issue May 7, 2022

type-flexible FFTs? JuliaMath/AbstractFFTs.jl#34

Open

milankl added the enhancement New feature or request label May 10, 2022

milankl mentioned this issue Jun 13, 2022

test with posits SpeedyWeather/SpeedyWeather.jl#77

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

FFTs with posits #64

FFTs with posits #64

milankl commented May 7, 2022

milankl commented May 10, 2022 •

edited

Loading

milankl commented Jun 10, 2022

milankl commented Jun 10, 2022

tkgunaratne commented Jun 10, 2022

giordano commented Jun 10, 2022

tkgunaratne commented Jun 10, 2022 •

edited

Loading

milankl commented Jun 11, 2022

milankl commented Jun 11, 2022

giordano commented Jun 11, 2022

milankl commented Jun 13, 2022

milankl commented Jun 13, 2022 •

edited

Loading

tkgunaratne commented Jun 13, 2022 •

edited

Loading

giordano commented Jun 13, 2022

tkgunaratne commented Jun 22, 2022

milankl commented Jun 23, 2022

milankl commented Jun 23, 2022

giordano commented Jun 23, 2022 •

edited

Loading

tkgunaratne commented Aug 20, 2022

milankl commented Aug 20, 2022

tkgunaratne commented Oct 11, 2022 via email

milankl commented Oct 11, 2022

tkgunaratne commented Oct 11, 2022 •

edited

Loading

milankl commented Oct 13, 2022

tkgunaratne commented Oct 13, 2022

tkgunaratne commented Apr 23, 2023

tkgunaratne commented Sep 17, 2024

FFTs with posits #64

FFTs with posits #64

Comments

milankl commented May 7, 2022

milankl commented May 10, 2022 • edited Loading

milankl commented Jun 10, 2022

milankl commented Jun 10, 2022

tkgunaratne commented Jun 10, 2022

giordano commented Jun 10, 2022

tkgunaratne commented Jun 10, 2022 • edited Loading

milankl commented Jun 11, 2022

milankl commented Jun 11, 2022

giordano commented Jun 11, 2022

milankl commented Jun 13, 2022

milankl commented Jun 13, 2022 • edited Loading

tkgunaratne commented Jun 13, 2022 • edited Loading

giordano commented Jun 13, 2022

tkgunaratne commented Jun 22, 2022

milankl commented Jun 23, 2022

milankl commented Jun 23, 2022

giordano commented Jun 23, 2022 • edited Loading

tkgunaratne commented Aug 20, 2022

milankl commented Aug 20, 2022

tkgunaratne commented Oct 11, 2022 via email

milankl commented Oct 11, 2022

tkgunaratne commented Oct 11, 2022 • edited Loading

milankl commented Oct 13, 2022

tkgunaratne commented Oct 13, 2022

tkgunaratne commented Apr 23, 2023

tkgunaratne commented Sep 17, 2024

milankl commented May 10, 2022 •

edited

Loading

tkgunaratne commented Jun 10, 2022 •

edited

Loading

milankl commented Jun 13, 2022 •

edited

Loading

tkgunaratne commented Jun 13, 2022 •

edited

Loading

giordano commented Jun 23, 2022 •

edited

Loading

tkgunaratne commented Oct 11, 2022 •

edited

Loading