1// Copyright ©2017 The Gonum Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package spectral
6
7import (
8	"math"
9
10	"gonum.org/v1/gonum/graph"
11	"gonum.org/v1/gonum/mat"
12)
13
14// Laplacian is a graph Laplacian matrix.
15type Laplacian struct {
16	// Matrix holds the Laplacian matrix.
17	mat.Matrix
18
19	// Nodes holds the input graph nodes.
20	Nodes []graph.Node
21
22	// Index is a mapping from the graph
23	// node IDs to row and column indices.
24	Index map[int64]int
25}
26
27// NewLaplacian returns a Laplacian matrix for the simple undirected graph g.
28// The Laplacian is defined as D-A where D is a diagonal matrix holding the
29// degree of each node and A is the graph adjacency matrix of the input graph.
30// If g contains self edges, NewLaplacian will panic.
31func NewLaplacian(g graph.Undirected) Laplacian {
32	nodes := graph.NodesOf(g.Nodes())
33	indexOf := make(map[int64]int, len(nodes))
34	for i, n := range nodes {
35		id := n.ID()
36		indexOf[id] = i
37	}
38
39	l := mat.NewSymDense(len(nodes), nil)
40	for j, u := range nodes {
41		uid := u.ID()
42		to := graph.NodesOf(g.From(uid))
43		l.SetSym(j, j, float64(len(to)))
44		for _, v := range to {
45			vid := v.ID()
46			if uid == vid {
47				panic("network: self edge in graph")
48			}
49			if uid < vid {
50				l.SetSym(indexOf[vid], j, -1)
51			}
52		}
53	}
54
55	return Laplacian{Matrix: l, Nodes: nodes, Index: indexOf}
56}
57
58// NewSymNormLaplacian returns a symmetric normalized Laplacian matrix for the
59// simple undirected graph g.
60// The normalized Laplacian is defined as I-D^(-1/2)AD^(-1/2) where D is a
61// diagonal matrix holding the degree of each node and A is the graph adjacency
62// matrix of the input graph.
63// If g contains self edges, NewSymNormLaplacian will panic.
64func NewSymNormLaplacian(g graph.Undirected) Laplacian {
65	nodes := graph.NodesOf(g.Nodes())
66	indexOf := make(map[int64]int, len(nodes))
67	for i, n := range nodes {
68		id := n.ID()
69		indexOf[id] = i
70	}
71
72	l := mat.NewSymDense(len(nodes), nil)
73	for j, u := range nodes {
74		uid := u.ID()
75		to := graph.NodesOf(g.From(uid))
76		if len(to) == 0 {
77			continue
78		}
79		l.SetSym(j, j, 1)
80		squdeg := math.Sqrt(float64(len(to)))
81		for _, v := range to {
82			vid := v.ID()
83			if uid == vid {
84				panic("network: self edge in graph")
85			}
86			if uid < vid {
87				to := g.From(vid)
88				k := to.Len()
89				if k < 0 {
90					k = len(graph.NodesOf(to))
91				}
92				l.SetSym(indexOf[vid], j, -1/(squdeg*math.Sqrt(float64(k))))
93			}
94		}
95	}
96
97	return Laplacian{Matrix: l, Nodes: nodes, Index: indexOf}
98}
99
100// NewRandomWalkLaplacian returns a damp-scaled random walk Laplacian matrix for
101// the simple graph g.
102// The random walk Laplacian is defined as I-D^(-1)A where D is a diagonal matrix
103// holding the degree of each node and A is the graph adjacency matrix of the input
104// graph.
105// If g contains self edges, NewRandomWalkLaplacian will panic.
106func NewRandomWalkLaplacian(g graph.Graph, damp float64) Laplacian {
107	nodes := graph.NodesOf(g.Nodes())
108	indexOf := make(map[int64]int, len(nodes))
109	for i, n := range nodes {
110		id := n.ID()
111		indexOf[id] = i
112	}
113
114	l := mat.NewDense(len(nodes), len(nodes), nil)
115	for j, u := range nodes {
116		uid := u.ID()
117		to := graph.NodesOf(g.From(uid))
118		if len(to) == 0 {
119			continue
120		}
121		l.Set(j, j, 1-damp)
122		rudeg := (damp - 1) / float64(len(to))
123		for _, v := range to {
124			vid := v.ID()
125			if uid == vid {
126				panic("network: self edge in graph")
127			}
128			l.Set(indexOf[vid], j, rudeg)
129		}
130	}
131
132	return Laplacian{Matrix: l, Nodes: nodes, Index: indexOf}
133}
134