IntegralEquations / HMatrices.jl

A Julia library for hierarchical matrices
MIT License
39 stars 3 forks source link

Creating a `KernelMatrix` with custom row & column elements #53

Closed tpgillam closed 5 months ago

tpgillam commented 5 months ago

In Kernel matrices documentation it's stated that elements can be of any type:

The row and columns elements may be as simple as points in $\mathbb{R}^d$ (as is the case for Nyström methods), but they can also be more complex objects such as triangles or basis functions –- the only thing required is that f(X[i],Y[j]) make sense.

However, I run into problems constructing ClusterTrees if I try to exercise this ability. Some assorted observations:

Am I approaching this case incorrectly? Or is this a bug? In any case, I'm looking forward to using HMatrices, as it should be directly applicable to my problem if it can handle this case. Thanks!

maltezfaria commented 5 months ago

This is certainly a bug. As you correctly pointed out, this package became independent of WavePropBase at some point, and I was not sufficiently careful to test the changes thoroughly. And the documentation is certainly outdated in the places you mentioned!

Let me take a look at these issues, and I will get back to you soon with more information.

maltezfaria commented 5 months ago

If you have a MWE of what you are trying to achieve, that would be quite helpful in understanding/fixing the problem! Thanks.

maltezfaria commented 5 months ago

I made some small changes here and added a simple example to kernelmatrix.md of how you may achieve something similar to what you asked. Let me know if this is not it :-)

tpgillam commented 5 months ago

If you have a MWE

Here's a very minimal MWE. (In my case the type actually represents a some geometry of finite extent, and the kernel is a function of a pair of elements, in case that were to change anything)

using HMatrices
using StaticArrays

struct MyType
    a::Float64
    b::Float64
end

HMatrices.center(t::MyType) = SVector(t.a, t.b)

X = [MyType(x, 2 * x) for x in 1:10]
Y = X

f(x::MyType, y::MyType) = (x.a + x.b + y.a + y.b) / 4

K = KernelMatrix(f, X, Y)
assemble_hmatrix(K)

which gives

ERROR: LoadError: MethodError: no method matching ClusterTree(::Vector{…}, ::HyperRectangle{…}, ::UnitRange{…}, ::Vector{…}, ::Vector{…}, ::Nothing, ::Nothing)

Closest candidates are:
  ClusterTree(::Array{SVector{N, T}, 1}, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) where {N, T}
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/clustertree.jl:25
  ClusterTree(::Any, ::Any; copy_elements, threads)
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/clustertree.jl:123
  ClusterTree(::Any; ...)
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/clustertree.jl:123

Stacktrace:
 [1] ClusterTree(elements::Vector{MyType}, splitter::CardinalitySplitter; copy_elements::Bool, threads::Bool)
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/clustertree.jl:138
 [2] ClusterTree(elements::Vector{MyType}, splitter::CardinalitySplitter)
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/clustertree.jl:123
 [3] assemble_hmatrix(K::KernelMatrix{…}; atol::Int64, rank::Int64, rtol::Float64, kwargs::@Kwargs{})
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/hmatrix.jl:324
 [4] assemble_hmatrix(K::KernelMatrix{typeof(f), Vector{MyType}, Vector{MyType}, Float64})
   @ HMatrices ~/.julia/packages/HMatrices/prQKn/src/hmatrix.jl:313
 [5] top-level scope
   @ ~/Documents/hmoo.jl:17
 [6] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [7] top-level scope
   @ REPL[3]:1
tpgillam commented 5 months ago

I made some small changes here and added a simple example to kernelmatrix.md of how you may achieve something similar to what you asked. Let me know if this is not it :-)

Thank you! I've tried running the example you give there (against the current release), and I get the same ClusterTree error as I quoted above.

tpgillam commented 5 months ago

Ahh sorry, I just saw the map(center, ...) that you added when moving the function. Won't test it right now but I suspect that should do the trick.

maltezfaria commented 5 months ago

Thanks for the MWE. I just tested it, and it runs on the dev branch. I will probably get this merged today, together with some other small changes I want to add, and tag a minor version. Let me know if that does not fix your issues.

For reference, it was indeed the commit simplifying the ClusterTree structure that broke the functionality you wanted!

maltezfaria commented 5 months ago

Fixed in #52