Sound use of multiple dispatch?

General Usage Performance
1

Is the code below a sound use of multiple dispatch?
Are there cleaner or more sound ways to the same?

using StaticArrays

const Point1D = SVector{1,Float64}
const Point2D = SVector{2,Float64}

" Generates a one-dimensional uniform mesh "
function genMesh(N;a=Point1D(0.),b=Point1D(1.))
    if (a.x>=b.x)
        ErrorException("Wrong input")
    else 
        mesh = [a + (i-1)*(b-a)/N for i=1:N+1]
    end 
end 

" Generates a two-dimensional uniform mesh "
function genMesh(N,a::Point2D,b::Point2D)
    if ((a.x>=b.x)||(a.y>=b.y))
        ErrorException("Wrong input")
    else         
        mesh = [Point2D(a.x + (i-1)*(b.x-a.x)/N, b.y + (j-1)*(b.y-a.y)/N) for i=1:N+1,j=1:N+1]
    end 
end

with output

# Generate a one-dimensional mesh with 5 intervals  
mesh = genMesh(5)

and

# Generate a two-dimensional mesh with 5 intervals  
mesh = genMesh(10,Point2D(0,0.),Point2D(1.,1.))

Thanks.

2

If you remove the semicolon in the argument list of the first method, e.g.

function genMesh(N,a=Point1D(0.),b=Point1D(1.))

then yes. The reason is that Julia does not dispatch on keyword arguments.

3 Likes
3

I have a hunch this method is also intended to annotate a::Point1D = Point1D(0.), b::Point1D = Point1D(1.). Default values don’t annotate arguments so the first method is genMesh(N::Any, a::Any, b::Any), which allows unsupported inputs to end up at that method instead of throwing a MethodError.

Also a.x and b.x seems like a mistake, SVectors don’t have fields like that, let alone expose them as API. This bothered me so I looked at the docs and it seems to have slipped my mind that the properties x, y, z, and w are implemented in StaticArrays.jl/src/MVector.jl for the first 4 elements of SVectors or MVectors. Also worth mentioning Swizzle.jl for an extension of this syntax.

2 Likes
4

If OP wants that syntax, use FieldVector subtypes, FieldVector is also defined in StaticArrays:

julia> struct Point2D{T} <: FieldVector{2,T}
           x::T
           y::T
       end

julia> p = Point2D(0.0, 1.0)
2-element Point2D{Float64} with indices SOneTo(2):
 0.0
 1.0

julia> p.x
0.0

julia> p.y
1.0

ps: Maybe a cleaner way to dispatch there is to use:

julia> genmesh(::Type{Point1D}; args...) = "1D mesh!"
genmesh (generic function with 1 method)

julia> genmesh(::Type{Point2D}; args...) = "2D mesh!"
genmesh (generic function with 2 methods)

julia> genmesh(Point1D)
"1D mesh!"

julia> genmesh(Point2D)
"2D mesh!"

although this is almost the same as defining two different function names.

4 Likes
5

Thank you for drawing my attention to default values not annotating argument. Much appreciated.

6

I am not sure how to distinguish between less clean and more clean here. Could you pls. elaborate?

In a next step I would like to use multiple dispatch to construct the graph Laplacian diffusion operator on the mesh. Is multiply dispatch on the eltype of elements in the mesh a sound approach here?

7

The Swizzle project is no longer actively developed it seems. The developer is working on Finch.jl now.

8

The point is: is you function genmesh always requiring from the user the definition of a and b points?

If the answer is yes, than you can dispatch on the type of those points:

julia> struct Point1D{T} <: FieldVector{1,T}
           x::T
       end

julia> struct Point2D{T} <: FieldVector{2,T}
           x::T
           y::T
       end

julia> genmesh(N::Int, a::Point1D, b::Point1D) = "1D mesh!"
genmesh (generic function with 1 method)

julia> genmesh(N::Int, a::Point2D, b::Point2D) = "2D mesh!"
genmesh (generic function with 2 methods)

julia> a = Point1D(0.0); b = Point1D(1.0);

julia> genmesh(100, a, b)
"1D mesh!"

julia> a = Point2D(0.0, 0.0); b = Point2D(1.0, 1.0);

julia> genmesh(100, a, b)
"2D mesh!"

but if, as in your example, you have default values for these a, and b parameters, and providing them explicitly is optional, there is no way for dispatch to know if the user wants a 1D or 2D mesh, if it can provide only the number of points. Thus, then, you need to provide the dimension somehow. One way is to define two functions, genmesh1D, genmesh2D, other way is to provide two methods but require as a parameter the type of point (dimension) you want, which is what I did there.

Cleaner or not is something that you have to choose from a user point of view, in this case.

3 Likes
9

This is the xyzw/rgba syntax package:
GitHub - serenity4/Swizzles.jl: Adding swizzles to your vectors
That note you are referencing belongs to this package:
GitHub - willow-ahrens/Swizzles.jl: AbstracterArrays

1 Like
10

I’m skeptical of functions that generate dense grids. Grids should be lazy, a 1d grid should just be a range, and a 2d grid should be 2 ranges, etc. I’ve never seen the need to collect them into dense arrays, they waste time and memory.

11

Aha!

Thus you suggest something like

const Point2D = SVector{2,Float64}

" Generates a one-dimensional uniform mesh between point 0 and 1"
function genMesh(N::Int)
    mesh = range(0,1,length=N+1) 
end 

" Generates a one-dimensional uniform mesh between point a and b"
function genMesh(N::Int,a::Real,b::Real)
    if (a>=b)
        ErrorException("Wrong input")
    else 
        mesh = range(a,b,length=N+1)
    end 
end 

" Generates a two-dimensional uniform mesh with same number of elements in both directions "
function genMesh(N::Int,a::Point2D,b::Point2D)
    if ((a.x>=b.x)||(a.y>=b.y))
        ErrorException("Wrong input")
    else         
        mesh = (range(a.x,b.x,length=N+1),range(a.y,b.y,length=N+1))
    end 
end

with sample output

mesh1d = genMesh(5,1.,5.)

1.0:0.8:5.0

and

mesh2d = genMesh(5,Point2D(1.,1.),Point2D(5.,6.))

(1.0:0.8:5.0, 1.0:1.0:6.0) .

Correct?

In a next step, I (obviously) would like to evaluate functions on the mesh. My first draft is

" Generates a vector on a one-dimensional uniform mesh by nodal evaluation a function"
function genVector(mesh::StepRangeLen,sourceFct::Function)  
    f = map(sourceFct,mesh)
    return f      
end

" Generates a vector on a two-dimensional uniform mesh by nodal evaluation a function"
function genVector(mesh::Tuple,sourceFct::Function)  
    f = map(sourceFct,Iterators.product(mesh2d[1],mesh2d[2]))
    f = reshape(f,:,1)
    return f      
end

Two questions emerge:

Q1: what is an appropriate way to annotate the type of the mesh (first input arg)?

Q2: what is an appropriate way to annotate the type of the function (second input arg)?

Thx!

12

You need to read about broadcasting: More Dots: Syntactic Loop Fusion in Julia
Single- and multi-dimensional Arrays · The Julia Language
:wink:

1 Like
13

True, correct. Your input is much appreciated.

" Generates a vector on a two-dimensional uniform mesh by nodal evaluation a function"
function genVector(mesh::Tuple,sourceFct::Function)  
    f = broadcast(sourceFct,mesh2d[1]',mesh2d[2])
    f = reshape(f,:,1)
    return f      
end

Two questions emerge:

Q1: what is an appropriate way to annotate the type of the mesh (first input arg)?

Q2: what is an appropriate way to annotate the type of the function (second input arg)?

Thx again.

14

Personally, I wouldn’t create the function genVector at all. Broadcasting is so basic in Julia, I’d just use it inline:

vec(sourceFct.(mesh[2], mesh[1]'))

The idiomatic syntax is using the dot, not explicitly calling the broadcast function.

Futhermore, your current code doesn’t create a vector but a matrix of size Nx1, vec will create a vector, while reusing the data.

15

Clear, thank you so much!

Next is the question of generating the matrix. Below is my first attempt.

Q: what is an appropriate way to annotate the type of the input argument called mesh?

" Generates the one-dimensional diffusion or graph Laplacian matrix " 
function genStiffMat(mesh::StepRangeLen)
    Np1 = length(mesh)
    N = Np1-1; h = norm(mesh[end] - mesh[1])/N; h2 = h*h; 
    e = ones(Np1); #..note that ones(N+1,1) does *not* work here 
    A = Tridiagonal(-e[2:end], 2*e, -e[2:end]); 
    A = (1/h2)*A;    
    return A     
end

" Generates the two-dimensional diffusion or graph Laplacian matrix " 
function genStiffMat(mesh::Tuple)
    Nxp1 = length(mesh[1]); Nyp1 = length(mesh[2]);   
    #..contribution in x-direction 
    A1Dxx = genStiffMat(mesh[1])
    A2Dxx = kron(A1Dxx,I(Nyp1))
    #..contribution in y-direction
    A1Dyy = genStiffMat(mesh[2])
    A2Dyy = kron(I(Nxp1),A1Dyy)
    #..sum of x-direction and y-direction contribution 
    A = A2Dxx + A2Dyy
    return A 
end
16

I would probably write the functions signatures like this:

function genStiffMat(r::AbstractRange)

and

function genStiffMat(r1::AbstractRange, r2::AbstractRange)

instead of using a Tuple. Also, StepRangeLen seems too specific, and excludes ranges like 1:5, for example.

I don’t really have much or any experience with Laplacian matrices, but it seems to me that the output of the code has a lot of redundancy, using lots of memory for storing almost no information. Is there no more compact way of expressing the needed operation?

17

Many thanks again. Much appreciated.

1/ Q1: Does it make sense here to keep the mesh stored as a tuple (n-tuble in n degrees for n=1,2 and n=3) and dispatch on the length of the tuple (or NTuple)?

2/ Q2: To some extend, probably yes, the code contains some redundancy. The matrices are stored as sparse, removing redundancy to some degree. Possibly Kronecker.jl can be used in the future. Not sure how LightGraphs.jl stored its graph Laplacian. Not sure how removing redundancy here affects the code further downstream.

18

Yes, it does make sense, but as a matter of taste I would prefer giving them as separate arguments. This is a style choice up to you. It is not uncommon to support both, though (see e.g. size, zeros, etc.), and that is quite simple:

Alternative 1:

# Dispatch on tuple length
foo(t::Tuple{AbstractRange}) = ...
foo(t::NTuple{2, AbstractRange}) = ...

# Define input with several args:
foo(t::AbstractRange...) = foo(t)

Alternative 2:

# Dispatch on number of arguments
foo(t::AbstractRange) = ...
foo(t1::AbstractRange, t2::AbstractRange) = ...

# Define for tuple
foo(t::NTuple{N, AbstractRange}) where {N} = foo(t...)

Then you can call the function either way, and you just pick which implementation you prefer.

2 Likes
19

Thanks! Again, much appreciated.