suffix array

通常の接尾辞配列(suffix array, SA)について書いておきます。

suffix array とは

短く言うと、全てのsuffixを列挙し、それらを辞書順にソートし、suffixの開始位置の順番を配列に格納したものです。例えば文字列"banana"に対するsuffix arrayは以下の通り{5,3,1,0,4,2}になります。

開始位置 suffix
5 a
3 ana
1 anana
0 banana
4 na
2 nana

suffix arrayが何かについてのより詳しい説明は省きます。

一部の界隈においては、suffix arrayを作る前に入力文字列の末端に終端文字を必ず付ける習慣があります。終端文字は'$'などで表され、他のあらゆる文字より辞書順で若いとされます。その界隈の人が書いたっぽい解説記事では、あらゆる説明の際に文字列の末端に$が付けられています。実際にはべつに$を付けなくても上述のようにsuffix arrayは作れますし使えます。とはいえ実際のところ$を付けることにしたほうが実装がシンプルになるので、この記事においても$を付けることにします。例えば文字列"banana$"に対するsuffix arrayは以下の通り{6,5,3,1,0,4,2}になります。

開始位置 suffix
6 $
5 a$
3 ana$
1 anana$
0 banana$
4 na$
2 nana$
実装がシンプルになることの例

例えばC++で文字列をvector<uint8_t>で持つとします。すると例えば"ana"は{97,110,97}で"anana"は{97,110,97,110,97}になります。"ana"は"anana"の完全なprefixなので、この2つを辞書順で比較するためにはout of rangeを防ぐ処理が必要です。終端文字があればそれが不要になります。

#include <iostream>
#include <vector>
// 終端文字なし
int main(void) {
	const std::vector<uint8_t>ana{ 97,110,97 };
	const std::vector<uint8_t>anana{ 97,110,97,110,97 };
	for (int i = 0;; ++i) {
		if (i == ana.size()) { std::cout << "ana < anana" << std::endl; break; }
		if (i == anana.size()) { std::cout << "ana > anana" << std::endl; break; }
		if (ana[i] < anana[i]) { std::cout << "ana < anana" << std::endl; break; }
		if (ana[i] > anana[i]) { std::cout << "ana > anana" << std::endl; break; }
	}
	return 0;
}

#include <iostream>
#include <vector>
// 終端文字あり
int main(void) {
	const std::vector<uint8_t>ana{ 97,110,97,0 };
	const std::vector<uint8_t>anana{ 97,110,97,110,97,0 };
	for (int i = 0;; ++i) {
		if (ana[i] < anana[i]) { std::cout << "ana < anana" << std::endl; break; }
		if (ana[i] > anana[i]) { std::cout << "ana > anana" << std::endl; break; }
	}
	return 0;
}
ナイーブな構築方法

ナイーブな構築方法では、suffixを実際に全列挙してソートします。時間計算量はO(N^{2}logN)です。

std::vector<int> BruteForceSuffixArrayConstruction(const std::vector<uint8_t>& input) {

	//終端文字は0で、終端文字以外は正の数だと仮定する。
	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] > 0);
	assert(input.back() == 0);

	//suffixを全列挙してソートする。
	std::vector<std::vector<uint8_t>>suffix;
	for (int i = 0; i < input.size(); ++i) {
		std::vector<uint8_t> s;
		for (int j = i; j < input.size(); ++j)s.push_back(input[j]);
		suffix.push_back(s);
	}
	sort(suffix.begin(), suffix.end());

	//ソートしたときの順番を記録して返す。ここでは順番をsuffixの長さから求めている。
	std::vector<int>suffix_array;
	for (int i = 0; i < suffix.size(); ++i) {
		suffix_array.push_back(input.size() - suffix[i].size());
	}
	return suffix_array;
}
Larsson-Sadakane法による構築

アルゴリズムの説明は省きます。時間計算量はO(N(logN)^{2})です。http://mlab.cb.k.u-tokyo.ac.jp/~moris/lecture/2009-winter-experiments/schedule.htmの Larsson_Sadakane.cpp のリンク先に完動するコードが載っています。このコードはLarge-Scale Genome Sequence Processingという本に書かれているコードが元になっていると思われます。

Large-Scale Genome Sequence Processing

Large-Scale Genome Sequence Processing

SA-IS法による構築

これもアルゴリズムの説明は省きます。時間計算量はO(N)です。元論文のAppendix1に完動するコードが載っています。
github.com
下記はこのAppendix1と同様なコードです。GitHubリポジトリのライセンスと同じくCC BY-SA 3.0です。

#include<iostream>
#include<vector>
#include<algorithm>
#include<cassert>

bool TGet(const std::vector<uint8_t>&t, const int index) {
	return (t[index / 8] & (uint8_t(1) << (index % 8))) ? true : false;
}
void TSet(std::vector<uint8_t>&t, const int index, const bool value) {
	t[index / 8] = value ?
		(t[index / 8] | (uint8_t(1) << (index % 8))) :
		(t[index / 8] & (~(uint8_t(1) << (index % 8))));
}
bool IsLMS(const std::vector<uint8_t>&t, const int index) {
	return (index > 0) && TGet(t, index) && (!TGet(t, index - 1));
}
template<typename T>void GetBuckets(
	const T* input,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	for (int i = 0; i <= K; ++i)bucket[i] = 0;
	for (int i = 0; i < n; ++i)++bucket[input[i]];
	int sum = 0;
	for (int i = 0; i <= K; ++i) {
		sum += bucket[i];
		bucket[i] = end ? sum : (sum - bucket[i]);
	}
}
template<typename T>void InducedSortArrayL(
	const T* input,
	const std::vector<uint8_t>&t,
	int* suffix_array,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	GetBuckets<T>(input, bucket, n, K, end);
	for (int i = 0; i < n; ++i) {
		const int j = suffix_array[i] - 1;
		if ((j >= 0) && (!TGet(t, j)))suffix_array[bucket[input[j]]++] = j;
	}
}
template<typename T>void InducedSortArrayS(
	const T* input,
	const std::vector<uint8_t>&t,
	int* suffix_array,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	GetBuckets<T>(input, bucket, n, K, end);
	for (int i = n - 1; i >= 0; --i) {
		const int j = suffix_array[i] - 1;
		if ((j >= 0) && TGet(t, j))suffix_array[--bucket[input[j]]] = j;
	}
}
template<typename T>void SA_IS(const T* input, int* suffix_array, const int n, const int K) {

	std::vector<uint8_t>t((n + 7) / 8, 0);
	TSet(t, n - 2, false);
	TSet(t, n - 1, true);
	for (int i = n - 3; i >= 0; --i) {
		bool sets = (input[i] < input[i + 1]) ||
			(input[i] == input[i + 1] && TGet(t, i + 1));
		TSet(t, i,
			(input[i] < input[i + 1]) ||
			(input[i] == input[i + 1] && TGet(t, i + 1))
		);
	}

	//stage 1
	std::vector<int>bucket(K + 1, 0);
	GetBuckets<T>(input, bucket, n, K, true);
	for (int i = 0; i < n; ++i)suffix_array[i] = -1;
	for (int i = 1; i < n; ++i)if (IsLMS(t, i))suffix_array[--bucket[input[i]]] = i;
	InducedSortArrayL<T>(input, t, suffix_array, bucket, n, K, false);
	InducedSortArrayS<T>(input, t, suffix_array, bucket, n, K, true);
	int n1 = 0;
	for (int i = 0; i < n; ++i) {
		if (IsLMS(t, suffix_array[i]))suffix_array[n1++] = suffix_array[i];
	}
	for (int i = n1; i < n; ++i)suffix_array[i] = -1;
	int name = 0, prev = -1;
	for (int i = 0; i < n1; ++i) {
		int pos = suffix_array[i];
		bool diff = false;
		for (int d = 0; d < n; ++d) {
			if ((prev == -1) ||
				(input[pos + d] != input[prev + d]) ||
				(TGet(t, pos + d) != TGet(t, prev + d))) {
				diff = true;
				break;
			}
			else if (d > 0 && (IsLMS(t, pos + d) || IsLMS(t, prev + d)))break;
		}
		if (diff) {
			name++;
			prev = pos;
		}
		pos = ((pos % 2) == 0) ? (pos / 2) : ((pos - 1) / 2);
		suffix_array[n1 + pos] = name - 1;
	}
	for (int i = n - 1, j = n - 1; i >= n1; --i) {
		if (suffix_array[i] >= 0)suffix_array[j--] = suffix_array[i];
	}

	//stage 2
	int* suffix_array1 = suffix_array;
	int* s1 = suffix_array + n - n1;
	if (name < n1) {
		SA_IS<int>(s1, suffix_array1, n1, name - 1);
	}
	else {
		for (int i = 0; i < n1; ++i)suffix_array1[s1[i]] = i;
	}

	//stage 3
	GetBuckets<T>(input, bucket, n, K, true);
	for (int i = 1, j = 0; i < n; ++i) {
		if (IsLMS(t, i))s1[j++] = i;
	}
	for (int i = 0; i < n1; ++i)suffix_array1[i] = s1[suffix_array1[i]];
	for (int i = n1; i < n; ++i)suffix_array[i] = -1;
	for (int i = n1 - 1; i >= 0; --i) {
		const int j = suffix_array[i];
		suffix_array[i] = -1;
		suffix_array[--bucket[input[j]]] = j;
	}
	InducedSortArrayL<T>(input, t, suffix_array, bucket, n, K, false);
	InducedSortArrayS<T>(input, t, suffix_array, bucket, n, K, true);
}
std::vector<int> SAISSuffixArrayConstruction(const std::vector<uint8_t>& input) {

	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
	assert(input.back() == 0);

	std::vector<int>suffix_array(input.size(), 0);
	SA_IS<uint8_t>(&input[0], &suffix_array[0], input.size(), 256);
	return suffix_array;
}
int main(void) {
	std::vector<uint8_t>banana{ 98,97,110,97,110,97,0 };
	const auto suffix_array = SAISSuffixArrayConstruction(banana);
	for (const int x : suffix_array) std::cout << x << std::endl;
	return 0;
}
使いみち

リファレンス文字列のSuffix Arrayがあれば、クエリ文字列との完全一致箇所をsuffix array上の二分探索で求められます。例えば以下のように書けます。以下のコードの返り値は一致箇所の開始点を全列挙したものです。これもCC BY-SA 3.0としておきます。

std::vector<int> SuffixArrayMatch(
	const std::vector<uint8_t>& ref,
	const std::vector<uint8_t>& query) {

	for (const auto x : ref)assert(x != 0);
	for (const auto x : query)assert(x != 0);

	auto ref_ = ref;
	ref_.push_back(0);
	const auto suffix_array = SAISSuffixArrayConstruction(ref_);

	int L = 0, R = suffix_array.size() - 1;
	while (L + 1 < R) {
		const int mid = (L + R) / 2;
		bool query_is_front = true;
		for (int i = 0; i < query.size(); ++i) {
			if (query[i] < ref_[suffix_array[mid] + i]) {
				query_is_front = true;
				break;
			}
			else if (query[i] > ref_[suffix_array[mid] + i]) {
				query_is_front = false;
				break;
			}
		}
		if (query_is_front)R = mid;
		else L = mid;
	}

	std::vector<int>match_start_pos;
	for (int pos = R; pos < suffix_array.size(); pos++) {
		if (suffix_array[pos] > ref_.size() - query.size())break;
		for (int i = 0; i < query.size(); ++i) {
			if (query[i] != ref_[suffix_array[pos] + i]) {
				goto LABEL;
			}
		}
		match_start_pos.push_back(suffix_array[pos]);
	}

LABEL:
	sort(match_start_pos.begin(), match_start_pos.end());
	return match_start_pos;
}
圧縮全文索引について

上述の通り、suffix arrayがあれば全文一致検索クエリを高速に処理できます。しかしsuffix arrayはふつうリファレンス文字列よりも大きい（上のコードだと4倍）ので、非常に長い文字列に対しては使えません。

一般に、全文検索を高速化するための追加データのことを全文索引と呼び、元の文字列よりも小さな全文索引のことを特に圧縮全文索引と呼びます。この記事で扱ったsuffix arrayは全文索引ですが圧縮全文索引ではありません。圧縮全文索引としてFM-index, (FM-index以外の)compressed suffix array(CSA), LZ圧縮ベースのindexingなどが知られています。FM-indexが特に高性能であると言われていてよく用いられます。
d.hatena.ne.jp
例えば上の記事でもそう書かれています。余談ですが、CSAという単語は(FM-index以外の)特定のアルゴリズムを示す場合と、FM-indexを含むある種のコンセプトを示す場合に分かれます。

特定のタスクにおいてCSAのほうが高性能だったと言っている論文もあって、「fm-index vs compressed suffix array」とかでググると見つかります。
Practical aspects of Compressed Suffix Arrays and FM-Index in Searching DNA Sequences - Semantic Scholar