Line data    Source code

       1             : """
       2             :     L = umbrella(mesh[; symmetric=true, typ=Float64])
       3             : 
       4             : Get uniform Laplacian (Leif Kobbelt's *umbrella operator*) either as
       5             : `symmetric` *sparse* matrix (i.e., with vertex [`degree`](@ref)s on
       6             : matrix diagonal) or with row weighting by inverse vertex degrees
       7             : (i.e., with unit-diagonal). The optional keyword argument `typ`
       8             : determines the `SparseMatrix{typ}`.
       9             : 
      10             : !!! note
      11             :     This implementation does *not* treat the boundary as a
      12             :     univariate polygon.
      13             : 
      14             : See also [`graphlaplacian`](@ref), [`adjacencymatrix`](@ref)
      15             : """
      16           4 : function umbrella(mesh::Mesh; symmetric::Bool=true, typ=Float64)
      17           2 :     A = adjacencymatrix(mesh; one=typ(1))
      18           2 :     V = Diagonal(vec(sum(A; dims=2)))
      19             : 
      20           2 :     U = A - V
      21             : 
      22           2 :     if symmetric
      23           1 :         U
      24             :     else
      25           1 :         V \ U
      26             :     end
      27             : end
      28             : 
      29             : abstract type GraphLaplacianEdgeType end
      30           4 : struct Undirected <: GraphLaplacianEdgeType end
      31           1 : struct InEdges <: GraphLaplacianEdgeType end
      32           1 : struct OutEdges <: GraphLaplacianEdgeType end
      33             : 
      34             : abstract type GraphLaplacianMode end
      35           4 : struct Plain <: GraphLaplacianMode end
      36           1 : struct Normalized <: GraphLaplacianMode end
      37           1 : struct NormalizedSymmetric <: GraphLaplacianMode end
      38             : 
      39             : """
      40             :     L = graphlaplacian(mesh[; edges = Undirected(),
      41             :                               mode = Plain(),
      42             :                               typ=Float64)
      43             : 
      44             : Get Graph Laplacian `L` with
      45             : 
      46             : - `edges = Undirected()`: ignore edge orientation,
      47             : - `edged ∈ [InEdges(), OutEdges()]`: consider only given orientation.
      48             : 
      49             : and
      50             : 
      51             : - `mode = Plain()` graph Laplacian ``D - A`` (for adjacency matrix ``A``
      52             :    and diagonal matrix with vertex [`degree`](@ref)s ``D``, equals
      53             :    `-umbrella(mesh; symmetric=true)`),
      54             : - `mode = Normalized()` normalized graph Laplacian ``D^{-1}\\,(D-A)``
      55             :   (equals `-umbrella(mesh; symmetric=false)`),
      56             : - `mode = NormalizedSymmetric()` symmetrically normalized graph Laplacian
      57             :   ``I\\,-\\,D^{-\\tfrac{1}{2}}\\,(D-A)\\,D^{-\\tfrac{1}{2}}``.
      58             : 
      59             : See also [`umbrella`](@ref), [`adjacencymatrix`](@ref)
      60             : """
      61          12 : function graphlaplacian(mesh::Mesh;
      62             :                         edges::E=Undirected(), mode::M=Plain(),
      63             :                         typ=Float64) where {E<:GraphLaplacianEdgeType,
      64             :                                             M<:GraphLaplacianMode}
      65           6 :     A, D = _graphlaplacian(mesh, edges, typ)
      66           6 :     _graphlaplacianmode(A, D, mode)
      67             : end
      68             : 
      69           4 : function _graphlaplacian(mesh::Mesh, ::Undirected, typ)
      70           4 :     A = adjacencymatrix(mesh; one=one(typ))
      71           4 :     D = Diagonal(vec(sum(A; dims=2)))
      72             : 
      73           4 :     A, D
      74             : end
      75             : 
      76           1 : function _graphlaplacian(mesh::Mesh, ::OutEdges, typ)
      77           1 :     A = adjacencymatrix(mesh; one=one(typ), boundary=false)
      78           1 :     D = Diagonal(vec(sum(A, dims=2)))
      79           1 :     A, D
      80             : end
      81             : 
      82           1 : function _graphlaplacian(mesh::Mesh, ::InEdges, typ)
      83           1 :     A = adjacencymatrix(mesh; one=one(typ), boundary=false)'
      84             : 
      85           1 :     D = Diagonal(vec(sum(A, dims=1)))
      86           1 :     A, D
      87             : end
      88             : 
      89           4 : _graphlaplacianmode(A, D, ::Plain) = D - A
      90             : 
      91           1 : _graphlaplacianmode(A, D, ::Normalized) = D \ (D - A)
      92             : 
      93           1 : function _graphlaplacianmode(A, D, ::NormalizedSymmetric)
      94           2 :     D2 = Diagonal(1 ./ sqrt.(diag(D)))
      95           1 :     I - D2 * A * D2
      96             : end