はじめに

「色々な機械学習処理をブラウザ上で試せるサイトを作った」中で実装したモデルの解説の六回目です。

今回はOPTICSの実装について解説します。

デモはこちらから。（TaskをClusteringにして、ModelのOPTICSを選択）
実際のコードはoptics.jsにあります。

なお、可視化部分については一切触れません。

概説

OPTICSについては、Wikipediaの擬似コードをそのままコードに起こしただけとなります。

ただし、内部で優先度付きキューを使用する必要があるため、その実装も行っています。

優先度付きキュー

優先度付きキューの実装はシンプルで、以下の関数を用意します。

push
データを追加する。引数は値と優先度の二つ。
move
データの優先度を変更する。
shift
優先度の高いデータを取得する。

注意すべき点としては、優先度に距離を直接設定しているため、「優先度の値が小さい方が優先度が高い」ということです。

内部実装は、値と優先度のペアの配列を保持するようにしています。

コード

優先度付きキュー

class PriorityQueue {
	constructor(arr) {
		this._value = arr || []
	}

	get length() {
		return this._value.length;
	}

	_sort() {
		this._value.sort((a, b) => a[1] - b[1])
	}

	push(value, priority) {
		this._value.push([value, priority])
		this._sort()
	}

	move(value, priority) {
		for (let i = 0; i < this.length; i++) {
			if (this._value[i][0] === value) {
				this._value[i][1] = priority;
				this._sort()
				return
			}
		}
		this.push(value, priority)
	}

	shift() {
		const [value, priority] = this._value.shift();
		return value;
	}
}

OPTICS

class OPTICS {
	// https://en.wikipedia.org/wiki/OPTICS_algorithm
	constructor(eps = Infinity, minPts = 5, metric = 'euclid') {
		this._eps = eps;
		this._minPts = minPts;

		this._metric = metric
		switch (this._metric) {
		case 'euclid':
			this._d = (a, b) => Math.sqrt(a.reduce((s, v, i) => s + (v - b[i]) ** 2, 0));
			break
		case 'manhattan':
			this._d = (a, b) => a.reduce((s, v, i) => s + Math.abs(v - b[i]), 0)
			break
		case 'chebyshev':
			this._d = (a, b) => Math.max(...a.map((v, i) => Math.abs(v - b[i])))
			break;
		}
	}

	fit(datas) {
		const n = datas.length;
		const d = Array(n);
		for (let i = 0; i < n; d[i++] = Array(n));
		for (let i = 0; i < n; i++) {
			for (let j = 0; j < i; j++) {
				const v = this._d(datas[i], datas[j]);
				d[i][j] = d[j][i] = v;
			}
		}
		const getNeighbors = (i) => {
			const neighbors = [];
			for (let k = 0; k < n; k++) {
				if (d[i][k] < this._eps) neighbors.push(k);
			}
			return neighbors
		}
		const coreDist = (i) => {
			const neighbors = getNeighbors(i).map(k => d[i][k]);
			if (neighbors.length <= this._minPts) return null;
			neighbors.sort((a, b) => a - b);
			return neighbors[this._minPts];
		}
		const reachabilityDist = (o, i) => {
			const cd = coreDist(i);
			if (cd === null) return cd;
			return Math.max(cd, d[i][o])
		}

		const processed = Array(n).fill(false);
		const rd = Array(n).fill(null);
		const update = (n, p, seeds) => {
			const cd = coreDist(p);
			if (cd === null) return
			for (const o of n) {
				if (processed[o]) continue
				const nrd = Math.max(cd, d[p][o]);
				if (rd[o] === null) {
					rd[o] = nrd;
					seeds.push(o, nrd);
				} else if (nrd < rd[o]) {
					rd[o] = nrd;
					seeds.move(o, nrd)
				}
			}
		}

		this._core_distance = [];
		for (let p = 0; p < n; p++) {
			if (processed[p]) continue;
			const neighbors = getNeighbors(p);
			processed[p] = true;
			const cd = coreDist(p)
			this._core_distance.push([p, cd])
			if (cd !== null) {
				const seeds = new PriorityQueue()
				update(neighbors, p, seeds)
				while (seeds.length > 0) {
					const q = seeds.shift();
					const nd = getNeighbors(q);
					processed[q] = true
					const cdq = coreDist(q)
					this._core_distance.push([q, cdq])
					if (cdq !== null) {
						update(nd, q, seeds);
					}
				}
			}
		}
	}

	predict(threshold = 0.1) {
		let c = 0;
		const n = this._core_distance.length
		const clusters = Array(n)
		for (let i = 0; i < n; i++) {
			const [k, d] = this._core_distance[i]
			clusters[k] = d > threshold ? ++c : c;
		}
		return clusters
	}
}