Started work on ExtractClusters algorithm for OPTICS proposed by Ankerst et al.

Author: mpraski

dbscan.go

@@ -95,8 +95,6 @@ func (c *dbscanClusterer) Learn(data [][]float64) error {
 	c.endWorkers()
 
 	c.v = nil
-	c.m = nil
-	c.w = nil
 	c.p = nil
 	c.r = nil
 
@@ -202,15 +200,15 @@ func (c *dbscanClusterer) nearest(p int, l *int, r *[]int) {
 
 	c.w.Add(c.s)
 
-	for i := 0; i < c.s; i++ {
-		if i == c.o {
+	for i := 0; i < c.l; i += c.f {
+		if c.l-i <= c.f {
 			b = c.l - 1
 		} else {
-			b = (i + 1) * c.f
+			b = i + c.f
 		}
 
 		c.j <- &rangeJob{
-			a: i * c.f,
+			a: i,
 			b: b,
 		}
 	}
@@ -233,6 +231,11 @@ func (c *dbscanClusterer) startWorkers() {
 
 func (c *dbscanClusterer) endWorkers() {
 	close(c.j)
+
+	c.j = nil
+
+	c.m = nil
+	c.w = nil
 }
 
 func (c *dbscanClusterer) nearestWorker() {
@@ -258,10 +261,8 @@ func (c *dbscanClusterer) numWorkers() int {
 		b = 10
 	} else if c.l < 100000 {
 		b = 100
-	} else if c.l < 1000000 {
-		b = 1000
 	} else {
-		b = 10000
+		b = 1000
 	}
 
 	if c.workers == 0 {

errors.go

@@ -10,4 +10,5 @@ var (
 	ErrZeroEpsilon    = errors.New("Epsilon cannot be 0")
 	ErrZeroMinpts     = errors.New("MinPts cannot be 0")
 	ErrZeroWorkers    = errors.New("Number of workers cannot be less than 0")
+	ErrZeroXi         = errors.New("Xi cannot be 0")
 )

optics.go

@@ -5,11 +5,19 @@ import (
 	"sync"
 )
 
+type steepDownArea struct {
+	start, end int
+	mib        float64
+}
+
 type opticsClusterer struct {
 	minpts  int
 	workers int
 	eps     float64
 
+	// 1 - xi
+	x float64
+
 	distance DistanceFunc
 
 	// slices holding the cluster mapping and sizes
@@ -38,7 +46,7 @@ type opticsClusterer struct {
 }
 
 /* Implementation of OPTICS algorithm with concurrent nearest neighbour computation */
-func OPTICS(minpts int, eps float64, workers int, distance DistanceFunc) (HardClusterer, error) {
+func OPTICS(minpts int, eps, xi float64, workers int, distance DistanceFunc) (HardClusterer, error) {
 	if minpts < 1 {
 		return nil, ErrZeroMinpts
 	}
@@ -51,6 +59,10 @@ func OPTICS(minpts int, eps float64, workers int, distance DistanceFunc) (HardCl
 		return nil, ErrZeroEpsilon
 	}
 
+	if xi <= 0 {
+		return nil, ErrZeroXi
+	}
+
 	var d DistanceFunc
 	{
 		if distance != nil {
@@ -64,6 +76,7 @@ func OPTICS(minpts int, eps float64, workers int, distance DistanceFunc) (HardCl
 		minpts:   minpts,
 		workers:  workers,
 		eps:      eps,
+		x:        1 - xi,
 		distance: d,
 	}, nil
 }
@@ -103,13 +116,14 @@ func (c *opticsClusterer) Learn(data [][]float64) error {
 	c.endWorkers()
 
 	c.v = nil
-	c.re = nil
-	c.so = nil
-	c.m = nil
-	c.w = nil
 	c.p = nil
 	c.r = nil
 
+	c.extract()
+
+	c.re = nil
+	c.so = nil
+
 	c.mu.Unlock()
 
 	return nil
@@ -230,6 +244,128 @@ func (c *opticsClusterer) update(p int, l *int, r *[]int, q *priorityQueue) {
 	}
 }
 
+func (c *opticsClusterer) extract() {
+	var (
+		i, e  int
+		mib   float64
+		areas []*steepDownArea = make([]*steepDownArea, 0)
+	)
+
+	for i < c.l-1 {
+		if c.re[c.so[i]] == nil || c.re[c.so[i+1]] == nil {
+			continue
+		}
+
+		mib = math.Max(mib, c.re[c.so[i]].p)
+
+		if c.isSteepDown(i, &e) {
+			as := areas[:0]
+			for j := 0; j < len(areas); j++ {
+				if c.re[c.so[areas[j].start]].p*c.x < mib {
+					continue
+				}
+
+				as = append(as, &steepDownArea{
+					start: areas[j].start,
+					end:   areas[j].end,
+					mib:   math.Max(areas[j].mib, mib),
+				})
+			}
+			areas = as
+
+			areas = append(areas, &steepDownArea{
+				start: i,
+				end:   e,
+			})
+
+			i = e + 1
+			mib = c.re[c.so[i]].p
+
+		} else if c.isSteepUp(i, &e) {
+			as := areas[:0]
+			for j := 0; j < len(areas); j++ {
+				if c.re[c.so[areas[j].start]].p*c.x < mib {
+					continue
+				}
+
+				as = append(as, &steepDownArea{
+					start: areas[j].start,
+					end:   areas[j].end,
+					mib:   math.Max(areas[j].mib, mib),
+				})
+			}
+			areas = as
+
+			i = e + 1
+			mib = c.re[c.so[i]].p
+
+			for j := 0; j < len(areas); j++ {
+				// check for cluster satisfying conditions outlined by Ankerst at al.
+			}
+		} else {
+			i++
+		}
+	}
+}
+
+func (c *opticsClusterer) isSteepDown(i int, e *int) bool {
+	if c.re[c.so[i]].p*c.x > c.re[c.so[i+1]].p {
+		return false
+	}
+
+	var counter, j int = 0, i + 1
+
+	*e = j
+
+	for {
+		if c.re[c.so[j]].p < c.re[c.so[j+1]].p {
+			break
+		}
+
+		if c.re[c.so[j]].p*c.x <= c.re[c.so[j+1]].p {
+			*e = j
+			counter = 0
+		} else {
+			counter++
+		}
+
+		if counter > c.minpts {
+			break
+		}
+	}
+
+	return *e != j
+}
+
+func (c *opticsClusterer) isSteepUp(i int, e *int) bool {
+	if c.re[c.so[i]].p > c.re[c.so[i+1]].p*c.x {
+		return false
+	}
+
+	var counter, j int = 0, i + 1
+
+	*e = j
+
+	for {
+		if c.re[c.so[j]].p > c.re[c.so[j+1]].p {
+			break
+		}
+
+		if c.re[c.so[j]].p <= c.re[c.so[j+1]].p*c.x {
+			*e = j
+			counter = 0
+		} else {
+			counter++
+		}
+
+		if counter > c.minpts {
+			break
+		}
+	}
+
+	return *e != j
+}
+
 /* Divide work among c.s workers, where c.s is determined
  * by the size of the data. This is based on an assumption that neighbour points of p
  * are located in relatively small subsection of the input data, so the dataset can be scanned
@@ -244,15 +380,15 @@ func (c *opticsClusterer) nearest(p int, l *int, r *[]int) {
 
 	c.w.Add(c.s)
 
-	for i := 0; i < c.s; i++ {
-		if i == c.o {
+	for i := 0; i < c.l; i += c.f {
+		if c.l-i <= c.f {
 			b = c.l - 1
 		} else {
-			b = (i + 1) * c.f
+			b = i + c.f
 		}
 
 		c.j <- &rangeJob{
-			a: i * c.f,
+			a: i,
 			b: b,
 		}
 	}
@@ -275,6 +411,11 @@ func (c *opticsClusterer) startWorkers() {
 
 func (c *opticsClusterer) endWorkers() {
 	close(c.j)
+
+	c.j = nil
+
+	c.m = nil
+	c.w = nil
 }
 
 func (c *opticsClusterer) nearestWorker() {
@@ -300,10 +441,8 @@ func (c *opticsClusterer) numWorkers() int {
 		b = 10
 	} else if c.l < 100000 {
 		b = 100
-	} else if c.l < 1000000 {
-		b = 1000
 	} else {
-		b = 10000
+		b = 1000
 	}
 
 	if c.workers == 0 {