Open Jilhg opened 5 months ago
This code is quite old, but in case it's helpful to you https://github.com/JuliaApproximation/SingularIntegralEquations.jl/blob/master/src/Operators/Convolution.jl
Thank you, that's indeed what I was looking for.
I tried to use the package, but already the ConvolutionHelmholtz.jl example did not work. Is the SingularIntegralEquations repository discontinued?
I managed to include the Convolution.jl file into my own project. However, I had to change spaces slightly, for instance rangespace(C::ConcreteConvolution{Laurent{D},BT}) where {D,BT} = space(C.G)
to rangespace(C::ConcreteConvolution{Laurent{D,R},BT}) where {D,R,BT} = space(C.G)
. Apparently at least Laurent{PeriodicSegment{Float64},Float64}
is not identified by Laurent{D} where D
.
Btw fyi, after testing the code I realized the following: On a domain $(0,b)$ the Convolution operator of a Laurent function works, but if I represent the function using Fourier I got a wrong result! (I checked the code and somehow it only uses the odd coefficients of the function, i.e. only the cosine part is involved). Also, if use a domain $(a,b)$ the result for a Laurent function is shifted depending on $a$ (I think this is because the code misses the extra phase).
Anyway, I think I will be able to implement what I need using Convolution.jl as a template. So thank you very much.
Just for curiosity: Is there a way of mapping a function on a domain to a canonical domain? What I have in mind is to implement the convolution operator on the canonical domain and then use this to generalize to other domains.
Hi @Jilhg, sorry for the slow reply! It isn't being actively maintained, but it wouldn't take much to get the code to work again. That code was quite experimental, so it's indeed possible that not all mapped domains are covered and all Fourier/Laurent variants.
I think the ApproxFun functions tocanonical
and fromcanonical
functions implement this mapping for you.
Hi @MikaelSlevinsky, no problem, you have been a great help already. Thank you again! I checked the two functions tocanonical
and fromcanonical
and from my understanding they return the identity function on the domain/space (maybe I am using them incorrectly).
But anyway, I don't need them anymore because I simply implemented the operator on any PeriodicSegment by hand (which is probably also much more efficient).
Btw, in case you decide to revive SingularIntegralEquations.jl, I would be happy to contribute my code: I defined the convolution operator for PeriodicSegment of any real interval for both Laurent and Fourier spaces; however, only for the usual line integral, not for the singular integral (not sure what that even does).
In case you don't want to revive SingularIntegralEquations.jl, but you or the community thinks having convolution operators implemented in ApproxFun would be a useful addition to the project in general, I could also try to make a pull request to implement a convolution operator in one of the maintained repositories (I guess ApproxFunFourier.jl would be natural candidate here). I have the code and it’s short and it works well, so I am quite confident that I could do that. However, just to warn you, it would be my first time contributing to such a big project. If you prefer not to include convolution operators or have doubts for any reason, I fully respect that. I just want to offer :)
What do you think?
Hi @Jilhg, did you post your code somewhere?
I'm currently in the need of a convolution operator myself and have tried to combine the code from @MikaelSlevinsky and your comments into some working code, unfortunately the type signatures are throwing me for a loop.
I've managed to implement a getindex(C::ConcreteConvolution{Fourier{D,R},T}, k::Integer, j::Integer) where {D <: PeriodicSegment, R, T}
method with the correct behavior but unfortunately when I try to apply it to a function defined in domainspace(C)
I get the following error
ConcreteConvolution : Fourier(【-1.0,1.0❫) → Fourier(【-1.0,1.0❫)
Error showing value of type ConcreteConvolution{ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}, Int64}:
ERROR: Override getindex(::ConcreteConvolution{ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}, Int64}, ::Integer, ::Integer)
Stacktrace:
[1] defaultgetindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::Int64, j::Int64)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:361
[2] getindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::Int64, j::Int64)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:350
[3] getindex(V::ApproxFunBase.SubOperator{Int64, ConcreteConvolution{…}, Tuple{…}, Tuple{…}, Tuple{…}}, k::Int64, j::Int64)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/SubOperator.jl:313
[4] default_BandedMatrix(S::ApproxFunBase.SubOperator{Int64, ConcreteConvolution{…}, Tuple{…}, Tuple{…}, Tuple{…}})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:970
[5] BandedMatrix
@ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:878 [inlined]
[6] AbstractArray
@ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:958 [inlined]
[7] defaultgetindex(::ConcreteConvolution{…}, ::Val{…}, ::UnitRange{…}, ::UnitRange{…})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:397
[8] defaultgetindex
@ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:378 [inlined]
[9] getindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::UnitRange{Int64}, j::UnitRange{Int64})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:350
[10] show(io::IOContext{Base.TTY}, mimetype::MIME{Symbol("text/plain")}, B::Operator; header::Bool)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/show.jl:99
[11] show(io::IOContext{Base.TTY}, mimetype::MIME{Symbol("text/plain")}, B::Operator)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/show.jl:89
I'm assuming this might be because I don't quite understand why Fourier
gets expanded into ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}
under the hood.
I have used the following preliminary code so far. I checked that it works for both Laurent and Fourier spaces. I am planning on cleaning it up a bit and creating a PR in the next couple of weeks. Any feedback is welcome :)
export Convolution
abstract type Convolution{S,T} <: Operator{T} end
struct ConcreteConvolution{S<:Space,T} <: Convolution{S,T}
G::Fun{S,T}
end
Convolution(G::Fun{S,T}) where {S,T} = ConcreteConvolution(G)
ApproxFunBase.domain(C::ConcreteConvolution) = ApproxFunBase.domain(C.G)
ApproxFunBase.domainspace(C::ConcreteConvolution) = ApproxFunBase.space(C.G)
ApproxFunBase.rangespace(C::ConcreteConvolution) = ApproxFunBase.UnsetSpace()
ApproxFunBase.bandwidths(C::ConcreteConvolution) = error("No range space attached to Convolution")
ApproxFunBase.getindex(C::ConcreteConvolution,k::Integer,j::Integer) = error("No range space attached to Convolution")
ApproxFunBase.rangespace(C::ConcreteConvolution{Laurent{D,R}}) where {D,R} = ApproxFunBase.space(C.G)
ApproxFunBase.rangespace(C::ConcreteConvolution{Fourier{D,R}}) where {D,R} = ApproxFunBase.space(C.G)
ApproxFunBase.bandwidths(C::ConcreteConvolution{Laurent{D,R}}) where {D,R} = 0,0
ApproxFunBase.bandwidths(C::ConcreteConvolution{Fourier{D,R}}) where {D,R} = 1,1
function ApproxFunBase.getindex(C::ConcreteConvolution{Laurent{D,R},T},k::Integer,j::Integer) where {D<:PeriodicSegment,R,T}
fourier_index::Integer = if isodd(k) div(k-1,2) else -div(k,2) end
if k == j && k ≤ ncoefficients(C.G)
return (exp(-2*pi*1im/arclength(domain(C.G))*fourier_index*first(domain(C.G)))*arclength(domain(C.G))*C.G.coefficients[k])::T
else
return zero(T)
end
end
function ApproxFunBase.getindex(C::ConcreteConvolution{Fourier{D,R},T},k::Integer,j::Integer) where {D<:PeriodicSegment,R,T}
fourier_index::Integer = if isodd(k) div(k-1,2) else div(k,2) end
if k == 1 && k ≤ ncoefficients(C.G)
if j==k
return (arclength(domain(C.G))*C.G.coefficients[1])::T
else
return zero(T)
end
elseif 2*fourier_index ≤ ncoefficients(C.G)
Gs::T = if 2*fourier_index ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index] else 0.0::T end # sine coefficient
Gc::T = if 2*fourier_index+1 ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index+1] else 0.0::T end # cosine coefficient
phase::T = 2*pi/arclength(domain(C.G))*fourier_index*first(domain(C.G))
if iseven(k) && j==k
return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T
elseif iseven(k) && j==k+1
return (arclength(domain(C.G))*(Gc*sin(phase)+Gs*cos(phase))/2)::T
elseif isodd(k) && j==k
return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T
elseif isodd(k) && j==k-1
return (arclength(domain(C.G))*(-Gc*sin(phase)-Gs*cos(phase))/2)::T
else
return zero(T)
end
else
return zero(T)
end
end
Just copy this code into any file MyConvolution.jl
and then put include("MyConvolution.jl")
in your code. The Convolution operator is then constructed for instance as
T = Convolution(Fun(sin,Fourier()))
which you can apply to a function normally.
Hi @Jilhg, thank you for the swift reply!
comparing to your code, it turns out I forgot to import ApproxFunBase.getindex
in my module!
Julia didn't find the implementation, prompting me to override getindex
(which I thought I did). I'll try out your code too!
Given a (periodic) function $f$ I would like to implement its convolution operator $C_f$, i.e. an operator that I can later apply on a function to give $C_f g = \int dy f(x-y)g(y)$.
After doing some search through the code, I don't think a convolution operator has been implemented in ApproxFun, though as a complete newbie, tbh, I am quite lost in the code. I could imagine a convolution was implemented as least implicitly, because I guess it is used for multiplying two functions over any Fourier space.
I managed to implement and explicit convolution for two functions on a Laurent() space. This is convenient, as the convolution basically just amounts to multiplying the coefficients of both functions:
It works, but it is of course not very flexible code (Note the extra phase has to be added in case the interval starts at $x_0 \neq 0$).
The reason I need to implement the convolution operator is that I would like to solve some linear equation that contains a convolution operator. So the above explicit implementation is not enough for me, I need to implement an operator that I can add and multiply to other operators.
From what I understand about the functionality of ApproxFun.jl I think a good implementation would look something like this:
I have the feeling all of these steps could be very easy to implement, but I am lacking the detailed understanding of how ApproxFun.jl works in the background. I would be great if someone could help me with that: How do I implement a custom operator in general? How do I take care of the conversion properly?
PS: In case others agree I would be great to add convolution operators to the general functionality of ApproxFun in the future. At least for me convolution operators appear frequently in scientific computing.