Stared work on gap statistic test for quality of k-means++ implementation

Author: mpraski

clusters.go

@@ -1,6 +1,8 @@
 package clusters
 
 import (
+	"math"
+
 	"gonum.org/v1/gonum/floats"
 )
 
@@ -23,10 +25,16 @@ type SCEvent struct {
 	Observation   []float64
 }
 
-/* Clusterer denotes a single operation of learning
+/* TestResult represents output of a test performed to measure quality of an algorithm. */
+type TestResult struct {
+	clusters, expected int
+}
+
+/* Clusterer denotes a operations of learning and testing
  * common for both Hard and Soft clusterers */
 type Clusterer interface {
-	Learn(data [][]float64) error
+	Learn([][]float64) error
+	Test([][]float64, ...interface{}) (*TestResult, error)
 }
 
 /* HardClusterer defines a set of operations for hard clustering algorithms */
@@ -80,4 +88,8 @@ var (
 	EuclideanDistance = func(a, b []float64) float64 {
 		return floats.Distance(a, b, 2)
 	}
+
+	EuclideanDistanceSquared = func(a, b []float64) float64 {
+		return math.Pow(floats.Distance(a, b, 2), 2)
+	}
 )

common.go

@@ -2,6 +2,10 @@ package clusters
 
 import (
 	"container/heap"
+	"math/rand"
+	"sync"
+
+	"gonum.org/v1/gonum/floats"
 )
 
 // struct denoting start and end indices of database portion to be scanned for nearest neighbours by workers in DBSCAN and OPTICS
@@ -63,3 +67,56 @@ func (pq *priorityQueue) Update(item *pItem, value int, priority float64) {
 	item.p = priority
 	heap.Fix(pq, item.i)
 }
+
+func wk(data [][]float64, centroids [][]float64, mapping []int) float64 {
+	var (
+		l  = float64(2 * len(data[0]))
+		wk = make([]float64, len(centroids))
+	)
+
+	for i := 0; i < len(mapping); i++ {
+		wk[mapping[i]-1] += EuclideanDistanceSquared(centroids[mapping[i]-1], data[i]) / l
+	}
+
+	return floats.Sum(wk)
+}
+
+func bounds(data [][]float64) []*[2]float64 {
+	var (
+		wg sync.WaitGroup
+
+		l = len(data[0])
+		r = make([]*[2]float64, l)
+	)
+
+	for i := 0; i < l; i++ {
+		r[i] = &[2]float64{
+			data[0][i],
+			data[0][i],
+		}
+	}
+
+	wg.Add(l)
+
+	for i := 0; i < l; i++ {
+		go func(n int) {
+			defer wg.Done()
+
+			for j := 0; j < len(data); j++ {
+				if data[j][n] < r[n][0] {
+					r[n][0] = data[j][n]
+				} else if data[j][n] > r[n][1] {
+					r[n][1] = data[j][n]
+				}
+			}
+		}(i)
+	}
+
+	wg.Wait()
+
+	return r
+}
+
+func uniform(data *[2]float64) float64 {
+	return rand.Float64()*(data[1]-data[0]) + data[0]
+}

dbscan.go

@@ -102,6 +102,14 @@ func (c *dbscanClusterer) Learn(data [][]float64) error {
 	return nil
 }
 
+func (c *dbscanClusterer) Test(data [][]float64, args ...interface{}) (*TestResult, error) {
+	if len(data) == 0 {
+		return nil, ErrEmptySet
+	}
+
+	return nil, nil
+}
+
 func (c *dbscanClusterer) Sizes() []int {
 	c.mu.RLock()
 	defer c.mu.RUnlock()

importer.go

@@ -8,14 +8,18 @@ import (
 	"strconv"
 )
 
-type Importer struct {
+type Importer interface {
+	Import(file string, start, end int) ([][]float64, error)
 }
 
-func NewImporter() *Importer {
-	return &Importer{}
+type csvImporter struct {
 }
 
-func (i *Importer) Import(file string, start, end, size int) ([][]float64, error) {
+func NewCsvImporter() Importer {
+	return &csvImporter{}
+}
+
+func (i *csvImporter) Import(file string, start, end int) ([][]float64, error) {
 	if start < 0 || end < 0 || start > end {
 		return [][]float64{}, ErrInvalidRange
 	}
@@ -28,7 +32,7 @@ func (i *Importer) Import(file string, start, end, size int) ([][]float64, error
 	defer f.Close()
 
 	var (
-		d = make([][]float64, 0, size)
+		d = make([][]float64, 0)
 		r = csv.NewReader(bufio.NewReader(f))
 		s = end - start + 1
 		g []float64

importer_test.go

@@ -10,11 +10,11 @@ const TOLERANCE = 0.000001
 func TestImportedLoadDataOfCorrectLengh(t *testing.T) {
 	var (
 		f = "data/test.csv"
-		i = NewImporter()
+		i = NewCsvImporter()
 		s = 3
 	)
 
-	d, e := i.Import(f, 0, 2, 3)
+	d, e := i.Import(f, 0, 2)
 	if e != nil {
 		t.Errorf("Error importing data: %s", e.Error())
 	}
@@ -27,15 +27,15 @@ func TestImportedLoadDataOfCorrectLengh(t *testing.T) {
 func TestImportedLoadCorrectData(t *testing.T) {
 	var (
 		f = "data/test.csv"
-		i = NewImporter()
+		i = NewCsvImporter()
 		s = [][]float64{
 			[]float64{0.1, 0.2, 0.3},
 			[]float64{0.4, 0.5, 0.6},
 			[]float64{0.7, 0.8, 0.9},
 		}
 	)
 
-	d, e := i.Import(f, 0, 2, 3)
+	d, e := i.Import(f, 0, 2)
 	if e != nil {
 		t.Errorf("Error importing data: %s", e.Error())
 	}
@@ -68,12 +68,12 @@ func fsliceEqual(a, b [][]float64) bool {
 func BenchmarkImport(b *testing.B) {
 	var (
 		f = "data/bus-stops.csv"
-		i = NewImporter()
+		i = NewCsvImporter()
 	)
 
 	b.ResetTimer()
 
-	_, e := i.Import(f, 4, 5, 15000)
+	_, e := i.Import(f, 4, 5)
 	if e != nil {
 		b.Errorf("Error importing data: %s", e.Error())
 	}

kmeans.go

@@ -1,6 +1,7 @@
 package clusters
 
 import (
+	"errors"
 	"math"
 	"math/rand"
 	"sync"
@@ -108,6 +109,52 @@ func (c *kmeansClusterer) Learn(data [][]float64) error {
 	return nil
 }
 
+func (c *kmeansClusterer) Test(data [][]float64, args ...interface{}) (*TestResult, error) {
+	if len(data) == 0 {
+		return nil, ErrEmptySet
+	}
+
+	clusters, ok := args[0].(int)
+	if !ok {
+		return nil, errors.New("Argument #0 is invalid")
+	}
+
+	var (
+		size   = len(data)
+		bounds = bounds(data)
+		wks    = make([]float64, clusters)
+		wkbs   = make([]float64, clusters)
+		sk     = make([]float64, clusters)
+		one    = make([]float64, clusters)
+		bwkbs  = make([]float64, clusters)
+	)
+
+	for i := 0; i < clusters; i++ {
+		c.Learn(data)
+
+		wks[i] = math.Log(wk(c.d, c.m, c.a))
+
+		for j := 0; j < clusters; j++ {
+			c.Learn(c.buildRandomizedSet(size, bounds))
+
+			bwkbs[j] = math.Log(wk(c.d, c.m, c.a))
+			one[j] = 1
+		}
+
+		wkbs[i] = floats.Sum(bwkbs) / float64(clusters)
+
+		floats.Scale(wkbs[i], one)
+		floats.Sub(bwkbs, one)
+		floats.Mul(bwkbs, bwkbs)
+
+		sk[i] = math.Sqrt(floats.Sum(bwkbs) / float64(clusters))
+	}
+
+	floats.Scale(math.Sqrt(1+(1/float64(clusters))), sk)
+
+	return nil, nil
+}
+
 func (c *kmeansClusterer) Sizes() []int {
 	c.mu.RLock()
 	defer c.mu.RUnlock()
@@ -324,3 +371,20 @@ func (c *kmeansClusterer) check() {
 
 	c.oldchanges = c.changes
 }
+
+func (c *kmeansClusterer) buildRandomizedSet(size int, bounds []*[2]float64) [][]float64 {
+	var (
+		l = len(bounds)
+		r = make([][]float64, size)
+	)
+
+	for i := 0; i < size; i++ {
+		r[i] = make([]float64, l)
+
+		for j := 0; j < l; j++ {
+			r[i][j] = uniform(bounds[j])
+		}
+	}
+
+	return r
+}

kmeans_test.go

@@ -1,6 +1,7 @@
 package clusters
 
 import (
+	"math"
 	"testing"
 )
 
@@ -11,10 +12,10 @@ func TestKmeansClusterNumerMatches(t *testing.T) {
 
 	var (
 		f = "data/bus-stops.csv"
-		i = NewImporter()
+		i = NewCsvImporter()
 	)
 
-	d, e := i.Import(f, 4, 5, 15000)
+	d, e := i.Import(f, 4, 5)
 	if e != nil {
 		t.Errorf("Error importing data: %s", e.Error())
 	}
@@ -32,3 +33,34 @@ func TestKmeansClusterNumerMatches(t *testing.T) {
 		t.Errorf("Number of clusters does not match: %d vs %d", len(c.Sizes()), C)
 	}
 }
+
+func TestKmeansUsingGapStatistic(t *testing.T) {
+	const (
+		C         = 8
+		THRESHOLD = 2
+	)
+
+	var (
+		f = "data/bus-stops.csv"
+		i = NewCsvImporter()
+	)
+
+	d, e := i.Import(f, 4, 5)
+	if e != nil {
+		t.Errorf("Error importing data: %s", e.Error())
+	}
+
+	c, e := KMeans(1000, C, EuclideanDistance)
+	if e != nil {
+		t.Errorf("Error initializing kmeans clusterer: %s", e.Error())
+	}
+
+	r, e := c.Test(d, 20)
+	if e != nil {
+		t.Errorf("Error running test: %s", e.Error())
+	}
+
+	if math.Abs(float64(r.clusters-r.expected)) > THRESHOLD {
+		t.Errorf("Discrepancy between numer of clusters and expectation: %d vs %d", r.clusters, r.expected)
+	}
+}

optics.go

@@ -135,6 +135,14 @@ func (c *opticsClusterer) Learn(data [][]float64) error {
 	return nil
 }
 
+func (c *opticsClusterer) Test(data [][]float64, args ...interface{}) (*TestResult, error) {
+	if len(data) == 0 {
+		return nil, ErrEmptySet
+	}
+
+	return nil, nil
+}
+
 func (c *opticsClusterer) Sizes() []int {
 	c.mu.RLock()
 	defer c.mu.RUnlock()