From 96d7b5d74ca6b5e68f2d5d7044d6c2f5e67a9fc9 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Sat, 28 Sep 2024 20:16:14 -0400 Subject: [PATCH] feat: cache plans for fft --- src/NeuralOperators.jl | 2 +- src/layers.jl | 23 ++++++++++++++------ src/transform.jl | 49 ++++++++++++++++++++++++++++++++++++------ 3 files changed, 61 insertions(+), 13 deletions(-) diff --git a/src/NeuralOperators.jl b/src/NeuralOperators.jl index 3c1ea89..869da88 100644 --- a/src/NeuralOperators.jl +++ b/src/NeuralOperators.jl @@ -3,7 +3,7 @@ module NeuralOperators using ArgCheck: @argcheck using ChainRulesCore: @non_differentiable using ConcreteStructs: @concrete -using FFTW: FFTW, irfft, rfft +using FFTW: FFTW, plan_rfft, plan_irfft using Random: Random, AbstractRNG using Static: StaticBool, False, True, known, static diff --git a/src/layers.jl b/src/layers.jl index dff7b13..b83307c 100644 --- a/src/layers.jl +++ b/src/layers.jl @@ -48,6 +48,14 @@ function LuxCore.initialparameters(rng::AbstractRNG, layer::OperatorConv) rng, eltype(layer.tform), out_chs, in_chs, layer.prod_modes)) end +function LuxCore.initialstates(::AbstractRNG, layer::OperatorConv) + fake_x = zeros(Float32, ntuple(Returns(1), ndims(layer.tform))..., 1) + plan_tform = plan_transform(layer.tform, fake_x, nothing) + x = transform(layer.tform, fake_x, plan_tform) + plan_inv_tform = plan_inverse(layer.tform, x, nothing, size(x)) + return (; plan_tform, plan_inv_tform) +end + function LuxCore.parameterlength(layer::OperatorConv) return layer.prod_modes * layer.in_chs * layer.out_chs end @@ -59,19 +67,21 @@ function OperatorConv( end function (conv::OperatorConv{True})(x::AbstractArray, ps, st) - return operator_conv(x, conv.tform, ps.weight), st + return operator_conv(x, conv.tform, ps.weight, st) end function (conv::OperatorConv{False})(x::AbstractArray, ps, st) N = ndims(conv.tform) xᵀ = permutedims(x, (ntuple(i -> i + 1, N)..., 1, N + 2)) - yᵀ = operator_conv(xᵀ, conv.tform, ps.weight) + yᵀ, stₙ = operator_conv(xᵀ, conv.tform, ps.weight, st) y = permutedims(yᵀ, (N + 1, 1:N..., N + 2)) - return y, st + return y, stₙ end -function operator_conv(x, tform::AbstractTransform, weights) - x_t = transform(tform, x) +function operator_conv(x, tform::AbstractTransform, weights, st) + plan_tform = plan_transform(tform, x, st.plan_tform) + x_t = transform(tform, x, plan_tform) + x_tr = truncate_modes(tform, x_t) x_p = apply_pattern(x_tr, weights) @@ -79,7 +89,8 @@ function operator_conv(x, tform::AbstractTransform, weights) x_padded = NNlib.pad_constant(x_p, expand_pad_dims(pad_dims), false; dims=ntuple(identity, ndims(x_p) - 2))::typeof(x_p) - return inverse(tform, x_padded, size(x)) + plan_inv_tform = plan_inverse(tform, x_padded, st.plan_inv_tform, size(x)) + return inverse(tform, x_padded, plan_inv_tform, size(x)), (; plan_tform, plan_inv_tform) end """ diff --git a/src/transform.jl b/src/transform.jl index 7fe53f2..10dd2ee 100644 --- a/src/transform.jl +++ b/src/transform.jl @@ -4,10 +4,22 @@ ## Interface - `Base.ndims(<:AbstractTransform)`: N dims of modes - - `transform(<:AbstractTransform, x::AbstractArray)`: Apply the transform to x - `truncate_modes(<:AbstractTransform, x_transformed::AbstractArray)`: Truncate modes that contribute to the noise - - `inverse(<:AbstractTransform, x_transformed::AbstractArray)`: Apply the inverse + +### Transform Interface + + - `plan_transform(<:AbstractTransform, x::AbstractArray, prev_plan)`: Construct a plan to + apply the transform to x. Might reuse the previous plan if possible + - `transform(<:AbstractTransform, x::AbstractArray, plan)`: Apply the transform to x using + the plan + +### Inverse Transform Interface + + - `plan_inverse(<:AbstractTransform, x_transformed::AbstractArray, prev_plan, M)`: + Construct a plan to apply the inverse transform to `x_transformed`. Might reuse the + previous plan if possible + - `inverse(<:AbstractTransform, x_transformed::AbstractArray, plan, M)`: Apply the inverse transform to `x_transformed` """ abstract type AbstractTransform{T} end @@ -22,7 +34,18 @@ end Base.ndims(T::FourierTransform) = length(T.modes) -transform(ft::FourierTransform, x::AbstractArray) = rfft(x, 1:ndims(ft)) +function plan_transform(ft::FourierTransform, x::AbstractArray, ::Nothing) + return plan_rfft(x, 1:ndims(ft)) +end + +function plan_transform(ft::FourierTransform, x::AbstractArray, prev_plan) + size(prev_plan) == size(x) && eltype(prev_plan) == eltype(x) && return prev_plan + return plan_transform(ft, x, nothing) +end + +@non_differentiable plan_transform(::Any...) + +transform(::FourierTransform, x::AbstractArray, plan) = plan * x function low_pass(ft::FourierTransform, x_fft::AbstractArray) return view(x_fft, map(d -> 1:d, ft.modes)..., :, :) @@ -30,7 +53,21 @@ end truncate_modes(ft::FourierTransform, x_fft::AbstractArray) = low_pass(ft, x_fft) -function inverse( - ft::FourierTransform, x_fft::AbstractArray{T, N}, M::NTuple{N, Int64}) where {T, N} - return real(irfft(x_fft, first(M), 1:ndims(ft))) +function plan_inverse(ft::FourierTransform, x_transformed::AbstractArray{T, N}, + ::Nothing, M::NTuple{N, Int64}) where {T, N} + return plan_irfft(x_transformed, first(M), 1:ndims(ft)) +end + +function plan_inverse(ft::FourierTransform, x_transformed::AbstractArray{T, N}, + prev_plan, M::NTuple{N, Int64}) where {T, N} + size(prev_plan) == size(x_transformed) && eltype(prev_plan) == eltype(x_transformed) && + return prev_plan + return plan_inverse(ft, x_transformed, nothing, M) +end + +@non_differentiable plan_inverse(::Any...) + +function inverse(::FourierTransform, x_transformed::AbstractArray{T, N}, plan, + ::NTuple{N, Int64}) where {T, N} + return real(plan * x_transformed) end