gonum.org/v1/gonum@v0.14.0/graph/network/diffusion_test.go (about)

     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  
     5  package network
     6  
     7  import (
     8  	"math"
     9  	"testing"
    10  
    11  	"gonum.org/v1/gonum/floats/scalar"
    12  	"gonum.org/v1/gonum/graph"
    13  	"gonum.org/v1/gonum/graph/simple"
    14  	"gonum.org/v1/gonum/graph/spectral"
    15  )
    16  
    17  var diffuseTests = []struct {
    18  	g []set
    19  	h map[int64]float64
    20  	t float64
    21  
    22  	wantTol float64
    23  	want    map[bool]map[int64]float64
    24  }{
    25  	{
    26  		g: grid(5),
    27  		h: map[int64]float64{0: 1},
    28  		t: 0.1,
    29  
    30  		wantTol: 1e-9,
    31  		want: map[bool]map[int64]float64{
    32  			false: {
    33  				A: 0.826684055, B: 0.078548060, C: 0.003858840, D: 0.000127487, E: 0.000003233,
    34  				F: 0.078548060, G: 0.007463308, H: 0.000366651, I: 0.000012113, J: 0.000000307,
    35  				K: 0.003858840, L: 0.000366651, M: 0.000018012, N: 0.000000595, O: 0.000000015,
    36  				P: 0.000127487, Q: 0.000012113, R: 0.000000595, S: 0.000000020, T: 0.000000000,
    37  				U: 0.000003233, V: 0.000000307, W: 0.000000015, X: 0.000000000, Y: 0.000000000,
    38  			},
    39  			true: {
    40  				A: 0.9063462486, B: 0.0369774705, C: 0.0006161414, D: 0.0000068453, E: 0.0000000699,
    41  				F: 0.0369774705, G: 0.0010670895, H: 0.0000148186, I: 0.0000001420, J: 0.0000000014,
    42  				K: 0.0006161414, L: 0.0000148186, M: 0.0000001852, N: 0.0000000016, O: 0.0000000000,
    43  				P: 0.0000068453, Q: 0.0000001420, R: 0.0000000016, S: 0.0000000000, T: 0.0000000000,
    44  				U: 0.0000000699, V: 0.0000000014, W: 0.0000000000, X: 0.0000000000, Y: 0.0000000000,
    45  			},
    46  		},
    47  	},
    48  	{
    49  		g: grid(5),
    50  		h: map[int64]float64{0: 1},
    51  		t: 1,
    52  
    53  		wantTol: 1e-9,
    54  		want: map[bool]map[int64]float64{
    55  			false: {
    56  				A: 0.2743435076, B: 0.1615920872, C: 0.0639346641, D: 0.0188054933, E: 0.0051023569,
    57  				F: 0.1615920872, G: 0.0951799548, H: 0.0376583937, I: 0.0110766934, J: 0.0030053582,
    58  				K: 0.0639346641, L: 0.0376583937, M: 0.0148997194, N: 0.0043825455, O: 0.0011890840,
    59  				P: 0.0188054933, Q: 0.0110766934, R: 0.0043825455, S: 0.0012890649, T: 0.0003497525,
    60  				U: 0.0051023569, V: 0.0030053582, W: 0.0011890840, X: 0.0003497525, Y: 0.0000948958,
    61  			},
    62  			true: {
    63  				A: 0.4323917545, B: 0.1660487336, C: 0.0270298904, D: 0.0029720194, E: 0.0003007247,
    64  				F: 0.1660487336, G: 0.0463974679, H: 0.0063556078, I: 0.0006056850, J: 0.0000589574,
    65  				K: 0.0270298904, L: 0.0063556078, M: 0.0007860810, N: 0.0000691647, O: 0.0000065586,
    66  				P: 0.0029720194, Q: 0.0006056850, R: 0.0000691647, S: 0.0000057466, T: 0.0000005475,
    67  				U: 0.0003007247, V: 0.0000589574, W: 0.0000065586, X: 0.0000005475, Y: 0.0000000555,
    68  			},
    69  		},
    70  	},
    71  	{
    72  		g: grid(5),
    73  		h: map[int64]float64{0: 1},
    74  		t: 10,
    75  
    76  		wantTol: 1e-9,
    77  		want: map[bool]map[int64]float64{
    78  			false: {
    79  				A: 0.0432375924, B: 0.0426071834, C: 0.0415872351, D: 0.0405673794, E: 0.0399371202,
    80  				F: 0.0426071834, G: 0.0419859658, H: 0.0409808885, I: 0.0399759024, J: 0.0393548325,
    81  				K: 0.0415872351, L: 0.0409808885, M: 0.0399998711, N: 0.0390189428, O: 0.0384127403,
    82  				P: 0.0405673794, Q: 0.0399759024, R: 0.0390189428, S: 0.0380620700, T: 0.0374707336,
    83  				U: 0.0399371202, V: 0.0393548325, W: 0.0384127403, X: 0.0374707336, Y: 0.0368885843,
    84  			},
    85  			true: {
    86  				A: 0.0532814862, B: 0.0594280160, C: 0.0462076361, D: 0.0330529557, E: 0.0211688130,
    87  				F: 0.0594280160, G: 0.0612529898, H: 0.0462850376, I: 0.0319891593, J: 0.0213123519,
    88  				K: 0.0462076361, L: 0.0462850376, M: 0.0340410963, N: 0.0229646704, O: 0.0152763556,
    89  				P: 0.0330529557, Q: 0.0319891593, R: 0.0229646704, S: 0.0153031853, T: 0.0103681461,
    90  				U: 0.0211688130, V: 0.0213123519, W: 0.0152763556, X: 0.0103681461, Y: 0.0068893147,
    91  			},
    92  		},
    93  	},
    94  	{
    95  		g: grid(5),
    96  		h: func() map[int64]float64 {
    97  			m := make(map[int64]float64, 25)
    98  			for i := int64(A); i <= Y; i++ {
    99  				m[i] = 1
   100  			}
   101  			return m
   102  		}(),
   103  		t: 1, // FIXME(kortschak): Low t used due to instability in mat.Exp.
   104  
   105  		wantTol: 1e-1, // FIXME(kortschak): High tolerance used due to instability in mat.Exp.
   106  		want: map[bool]map[int64]float64{
   107  			false: {
   108  				A: 1, B: 1, C: 1, D: 1, E: 1,
   109  				F: 1, G: 1, H: 1, I: 1, J: 1,
   110  				K: 1, L: 1, M: 1, N: 1, O: 1,
   111  				P: 1, Q: 1, R: 1, S: 1, T: 1,
   112  				U: 1, V: 1, W: 1, X: 1, Y: 1,
   113  			},
   114  			true: {
   115  				// Output from the python implementation associated with doi:10.1371/journal.pcbi.1005598.
   116  				A: 0.98264450473308107, B: 1.002568278028513, C: 0.9958911385307706, D: 1.002568278028513, E: 0.98264450473308107,
   117  				F: 1.002568278028513, G: 1.0075291695232433, H: 1.0038067383118021, I: 1.0075291695232433, J: 1.002568278028513,
   118  				K: 0.9958911385307706, L: 1.0038067383118021, M: 1.0001850837547184, N: 1.0038067383118021, O: 0.9958911385307706,
   119  				P: 1.002568278028513, Q: 1.0075291695232433, R: 1.0038067383118021, S: 1.0075291695232433, T: 1.002568278028513,
   120  				U: 0.98264450473308107, V: 1.002568278028513, W: 0.9958911385307706, X: 1.002568278028513, Y: 0.98264450473308107,
   121  			},
   122  		},
   123  	},
   124  	{
   125  		g: []set{
   126  			A: linksTo(B, C),
   127  			B: linksTo(D),
   128  			C: nil,
   129  			D: nil,
   130  			E: linksTo(F),
   131  			F: nil,
   132  		},
   133  		h: map[int64]float64{A: 1, E: 10},
   134  		t: 0.1,
   135  
   136  		wantTol: 1e-9,
   137  		want: map[bool]map[int64]float64{
   138  			false: {
   139  				A: 0.8270754166, B: 0.0822899600, C: 0.0863904410, D: 0.0042441824, E: 9.0936537654, F: 0.9063462346,
   140  			},
   141  			true: {
   142  				A: 0.9082331512, B: 0.0453361743, C: 0.0640616812, D: 0.0016012085, E: 9.0936537654, F: 0.9063462346,
   143  			},
   144  		},
   145  	},
   146  }
   147  
   148  func TestDiffuse(t *testing.T) {
   149  	for i, test := range diffuseTests {
   150  		g := simple.NewUndirectedGraph()
   151  		for u, e := range test.g {
   152  			// Add nodes that are not defined by an edge.
   153  			if g.Node(int64(u)) == nil {
   154  				g.AddNode(simple.Node(u))
   155  			}
   156  			for v := range e {
   157  				g.SetEdge(simple.Edge{F: simple.Node(u), T: simple.Node(v)})
   158  			}
   159  		}
   160  		for j, lfn := range []func(g graph.Undirected) spectral.Laplacian{spectral.NewLaplacian, spectral.NewSymNormLaplacian} {
   161  			normalize := j == 1
   162  			var wantTemp float64
   163  			for _, v := range test.h {
   164  				wantTemp += v
   165  			}
   166  			got := Diffuse(nil, test.h, lfn(g), test.t)
   167  			prec := 1 - int(math.Log10(test.wantTol))
   168  			for n := range test.g {
   169  				if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[normalize][int64(n)], test.wantTol, test.wantTol) {
   170  					t.Errorf("unexpected Diffuse result for test %d with normalize=%t:\ngot: %v\nwant:%v",
   171  						i, normalize, orderedFloats(got, prec), orderedFloats(test.want[normalize], prec))
   172  					break
   173  				}
   174  			}
   175  
   176  			if j == 1 {
   177  				continue
   178  			}
   179  
   180  			var gotTemp float64
   181  			for _, v := range got {
   182  				gotTemp += v
   183  			}
   184  			gotTemp /= float64(len(got))
   185  			wantTemp /= float64(len(got))
   186  			if !scalar.EqualWithinAbsOrRel(gotTemp, wantTemp, test.wantTol, test.wantTol) {
   187  				t.Errorf("unexpected total heat for test %d with normalize=%t: got:%v want:%v",
   188  					i, normalize, gotTemp, wantTemp)
   189  			}
   190  		}
   191  	}
   192  }
   193  
   194  var diffuseToEquilibriumTests = []struct {
   195  	g       []set
   196  	builder builder
   197  	h       map[int64]float64
   198  	damp    float64
   199  	tol     float64
   200  	iter    int
   201  
   202  	want   map[int64]float64
   203  	wantOK bool
   204  }{
   205  	{
   206  		g:       grid(5),
   207  		builder: simple.NewUndirectedGraph(),
   208  		h:       map[int64]float64{0: 1},
   209  		damp:    0.85,
   210  		tol:     1e-6,
   211  		iter:    1e4,
   212  
   213  		want: map[int64]float64{
   214  			A: 0.025000, B: 0.037500, C: 0.037500, D: 0.037500, E: 0.025000,
   215  			F: 0.037500, G: 0.050000, H: 0.050000, I: 0.050000, J: 0.037500,
   216  			K: 0.037500, L: 0.050000, M: 0.050000, N: 0.050000, O: 0.037500,
   217  			P: 0.037500, Q: 0.050000, R: 0.050000, S: 0.050000, T: 0.037500,
   218  			U: 0.025000, V: 0.037500, W: 0.037500, X: 0.037500, Y: 0.025000,
   219  		},
   220  		wantOK: true,
   221  	},
   222  	{
   223  		// Example graph from http://en.wikipedia.org/wiki/File:PageRanks-Example.svg 16:17, 8 July 2009
   224  		g: []set{
   225  			A: nil,
   226  			B: linksTo(C),
   227  			C: linksTo(B),
   228  			D: linksTo(A, B),
   229  			E: linksTo(D, B, F),
   230  			F: linksTo(B, E),
   231  			G: linksTo(B, E),
   232  			H: linksTo(B, E),
   233  			I: linksTo(B, E),
   234  			J: linksTo(E),
   235  			K: linksTo(E),
   236  		},
   237  		builder: simple.NewDirectedGraph(),
   238  		h: map[int64]float64{
   239  			A: 1. / 11.,
   240  			B: 1. / 11.,
   241  			C: 1. / 11.,
   242  			D: 1. / 11.,
   243  			E: 1. / 11.,
   244  			F: 1. / 11.,
   245  			G: 1. / 11.,
   246  			H: 1. / 11.,
   247  			I: 1. / 11.,
   248  			J: 1. / 11.,
   249  			K: 1. / 11.,
   250  		},
   251  		damp: 0.85,
   252  		tol:  1e-6,
   253  		iter: 1e4,
   254  
   255  		// This does not look like Page Rank because we do not
   256  		// do the random node hops. An alternative Laplacian
   257  		// value that does do that would replicate PageRank. This
   258  		// is left as an exercise for the reader.
   259  		want: map[int64]float64{
   260  			A: 0.227273,
   261  			B: 0.386364,
   262  			C: 0.386364,
   263  			D: 0.000000,
   264  			E: 0.000000,
   265  			F: 0.000000,
   266  			G: 0.000000,
   267  			H: 0.000000,
   268  			I: 0.000000,
   269  			J: 0.000000,
   270  			K: 0.000000,
   271  		},
   272  		wantOK: true,
   273  	},
   274  	{
   275  		g: []set{
   276  			A: linksTo(B, C),
   277  			B: linksTo(D, C),
   278  			C: nil,
   279  			D: nil,
   280  			E: linksTo(F),
   281  			F: nil,
   282  		},
   283  		builder: simple.NewDirectedGraph(),
   284  		h:       map[int64]float64{A: 1, E: -10},
   285  		tol:     1e-6,
   286  		iter:    3,
   287  
   288  		want: map[int64]float64{
   289  			A: 0, B: 0, C: 0.75, D: 0.25, E: 0, F: -10,
   290  		},
   291  		wantOK: true,
   292  	},
   293  	{
   294  		g: []set{
   295  			A: linksTo(B, C),
   296  			B: linksTo(D, C),
   297  			C: nil,
   298  			D: nil,
   299  			E: linksTo(F),
   300  			F: nil,
   301  		},
   302  		builder: simple.NewUndirectedGraph(),
   303  		h:       map[int64]float64{A: 1, E: -10},
   304  		damp:    0.85,
   305  		tol:     1e-6,
   306  		iter:    1e4,
   307  
   308  		want: map[int64]float64{
   309  			A: 0.25, B: 0.375, C: 0.25, D: 0.125, E: -5, F: -5,
   310  		},
   311  		wantOK: true,
   312  	},
   313  	{
   314  		g: []set{
   315  			A: linksTo(B),
   316  			B: linksTo(C),
   317  			C: nil,
   318  		},
   319  		builder: simple.NewUndirectedGraph(),
   320  		h:       map[int64]float64{B: 1},
   321  		iter:    1,
   322  		tol:     1e-6,
   323  		want: map[int64]float64{
   324  			A: 0.5, B: 0, C: 0.5,
   325  		},
   326  		wantOK: false,
   327  	},
   328  	{
   329  		g: []set{
   330  			A: linksTo(B),
   331  			B: linksTo(C),
   332  			C: nil,
   333  		},
   334  		builder: simple.NewUndirectedGraph(),
   335  		h:       map[int64]float64{B: 1},
   336  		iter:    2,
   337  		tol:     1e-6,
   338  		want: map[int64]float64{
   339  			A: 0, B: 1, C: 0,
   340  		},
   341  		wantOK: false,
   342  	},
   343  }
   344  
   345  func TestDiffuseToEquilibrium(t *testing.T) {
   346  	for i, test := range diffuseToEquilibriumTests {
   347  		g := test.builder
   348  		for u, e := range test.g {
   349  			// Add nodes that are not defined by an edge.
   350  			if g.Node(int64(u)) == nil {
   351  				g.AddNode(simple.Node(u))
   352  			}
   353  			for v := range e {
   354  				g.SetEdge(simple.Edge{F: simple.Node(u), T: simple.Node(v)})
   355  			}
   356  		}
   357  		var wantTemp float64
   358  		for _, v := range test.h {
   359  			wantTemp += v
   360  		}
   361  		got, ok := DiffuseToEquilibrium(nil, test.h, spectral.NewRandomWalkLaplacian(g, test.damp), test.tol*test.tol, test.iter)
   362  		if ok != test.wantOK {
   363  			t.Errorf("unexpected success value for test %d: got:%t want:%t", i, ok, test.wantOK)
   364  		}
   365  		prec := -int(math.Log10(test.tol))
   366  		for n := range test.g {
   367  			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.tol, test.tol) {
   368  				t.Errorf("unexpected DiffuseToEquilibrium result for test %d:\ngot: %v\nwant:%v",
   369  					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
   370  				break
   371  			}
   372  		}
   373  
   374  		var gotTemp float64
   375  		for _, v := range got {
   376  			gotTemp += v
   377  		}
   378  		gotTemp /= float64(len(got))
   379  		wantTemp /= float64(len(got))
   380  		if !scalar.EqualWithinAbsOrRel(gotTemp, wantTemp, test.tol, test.tol) {
   381  			t.Errorf("unexpected total heat for test %d: got:%v want:%v",
   382  				i, gotTemp, wantTemp)
   383  		}
   384  	}
   385  }
   386  
   387  type builder interface {
   388  	graph.Graph
   389  	graph.Builder
   390  }
   391  
   392  func grid(d int) []set {
   393  	dim := int64(d)
   394  	s := make([]set, dim*dim)
   395  	for i := range s {
   396  		s[i] = make(set)
   397  	}
   398  	for i := int64(0); i < dim*dim; i++ {
   399  		if i%dim != 0 {
   400  			s[i][i-1] = struct{}{}
   401  		}
   402  		if i/dim != 0 {
   403  			s[i][i-dim] = struct{}{}
   404  		}
   405  	}
   406  	return s
   407  }