【c/c++】Suffix ArrayとBWTに関するプログラムを作ってみた【c++17】

c/c++ プログラミング 日記

2018/6/25：concurrency::concurrent_vectorを使うべきだったがstd::vectorを使用しており、長い文字列の処理が不安定であったため、これを踏まえて全体に書き直し。


最近c++を触る機会ができたり、c++17の機能を使ってみたかったりしたので、課題にかこつけて取り組んでみました。

何を作ったか

Suffix Arrayと、Suffix ArrayからBWTの作成、BWTから元の文字列の復元をするプログラムを作りました。

Suffix Arrayとは

shogo82148.github.io

Burrows-Wheeler Transform(BWT)とは

d.hatena.ne.jp

プログラム

ソースコード

プロジェクト全体はGitHubに上げています。
github.com

#include <algorithm> //sort用
#include <execution> //algorithmの並列実行ポリシー
#include <ppl.h> //parallel_for用
#include <concurrent_vector.h> //並列

#include <iostream>
#include <string>
#include <tuple>

//単語を入れたらSuffixのベクターを返す
concurrency::concurrent_vector<std::tuple<int, std::string>>
MakeSuffix(const std::string &str
){
    using namespace std;

    concurrency::concurrent_vector<std::tuple<int, std::string>> v;
    //parallel_forで並列処理
    concurrency::parallel_for(size_t(0), str.length(), [&str, &v](size_t i){
        v.push_back(make_tuple(i, str.substr(i, str.length())));
    });

    return v;
}

//Suffixをソート
void
DictionaryOrderSort(concurrency::concurrent_vector<std::tuple<int, std::string>> &Suffix
){
    using namespace std;
    sort(execution::par_unseq,
         Suffix.begin(), Suffix.end(),
         [](tuple<int, string> &t1, tuple<int, string> &t2) { return (get<1>(t1) <= get<1>(t2)); }
        );
}

//BWT系列を作成
std::string
MakeBWT(const std::string OriginalString,
        const concurrency::concurrent_vector<std::tuple<int, std::string>> &SuffixArray
){
    using namespace std;

    string s(OriginalString.length(), '0'); //容量固定のstringを配列っぽく使用
    concurrency::parallel_for(size_t(0), OriginalString.length(), [&](size_t i){
        if(get<0>(SuffixArray[i]) == 0) s[i] = OriginalString[OriginalString.length()-1];
        else s[i] = OriginalString[get<0>(SuffixArray[i])-1];
    });

    return s;
}

std::string
ReconstructionFromBWT(const std::string &BWT,
                      const unsigned int limit = INT_MAX // 何文字目まで再生するかの指定
){
    using namespace std;
    //タグ付け
    vector<tuple<char, int>> v;
    string temp = BWT;
    size_t val; //BTWの行末をまず入れる
    for(size_t i = 0; i < temp.length(); i++){
        v.emplace_back(temp[i], i); //emplace_backの中に書いた要素でデフォルトコンストラクタが呼ばれてるそうな
        if(temp[i] == '$') val = i;
    }
    //安定ソートで並べ替え、法則を得る
    stable_sort(execution::par_unseq, v.begin(), v.end());

    //リミットまで復元
    const int lim = max((int)temp.length()-(int)limit-1, 0);
    int j;
    for(int i = (int)temp.length() -1; i >= lim; i--){
        temp[i] = BWT[val];
        j = 0;
        while(val != get<1>(v[j])) j++;
        val = j;
    }

    return temp.substr(max(lim, 0), temp.length());
}

int main() {
    using namespace std;

    //string str = "abca$";
    string str = "internationalization$"; //入力
    cout << "Suffix\n";
    for(auto t : v){
        cout << get<0>(t) << "\t:" << get<1>(t) << '\n';
    }

    DictionaryOrderSort(v); //ソート
    cout << "\nSuffix Array\n";
    for(auto t : v){
        cout << get<0>(t) << "\t:" << get<1>(t) << '\n';
    }

    string BWT = MakeBWT(str, v); //BWT取得
    cout << "\nGet BWT" << '\n';
    cout << BWT << '\n';

    cout << "\nDecode" << '\n';
    if(str == ReconstructionFromBWT(BWT)) cout << "success!";
    else cout << "failed";


    cout << flush;

    return 0;
}

出力

「internationalization」という単語についての出力です。

解説

重要なのは大きく以下の2点です。

  1. Suffix ArrayからBWT系列を求める場合の行末の扱い
  2. BWT系列から文字列を復元する方法
Suffix ArrayからBWT系列を求める場合の行末の扱い

Suffix ArrayからBWTを求める場合、行末を区別しなければソート順が乱れるパターンが有るため、何らかの記号を補っておく必要があります。単語の末尾に「$」を加えているのはそのためです。

BWT系列から文字列を復元する方

原理の説明は飛ばします。
BWT系列からの復元は、BWT系列の各文字に数字でタグを付け、

vector<tuple<char, int>> v;
string temp = get<1>(BWT);
for(size_t i = 0; i < temp.length(); i++){
    v.emplace_back(make_tuple(temp[i], i));
}

BWT文字列を辞書順に安定ソートし、

stable_sort(execution::par_unseq, v.begin(), v.end());

BWTキーから順に『キーと同値のタグを見つける→文字列の先頭にBWT文字列[タグ]を加える→見つけた場所の番号のタグが入っている場所をキーとする』という操作を繰り返していくことで、

int val = get<0>(BWT);
int j;
for(int i = (int)temp.length() -1; i >= lim; i--){
    temp[i] = get<1>(BWT)[val];
    j = 0;
    while(val != get<1>(v[j])) j++;
    val = j;
}

行うことができます。ただし復元は末尾から行われます。
今回のプログラムではMakeBWT関数で行末の位置を保持して返すようにしていますが、行末記号は文字列中に有るため、保持が無くても復元は行なえます。

工夫した点

parallel_forを使う

ppl.hのparallel_forを使って並列処理でSuffixを作りました。c++のラムダ式を触ってみたかったのですが、これで理解できたかなと思います。
場合によってはalgolithmのfor_eachにExecutionPolicyを指定する形とした方がppl.hに依存しなくなるので良いかなと思いました。……本音を言うとparallel_forの形が規格として取り込まれてくれるといいなと思います。
最初はparallel_forを使う時にconcurrent_vectorを使っていなかったせいで不安定になっていたのが反省点です。
参考資料：
parallel_forの書き方：
How to: Write a parallel_for Loop | Microsoft Docs
ラムダ式の書き方：
cpprefjp.github.io

std::sortの並列実行

c++17でalgorithmに並列実行オプションが追加されました。execution::par_unseqがそれです。c++17の機能、これしか使えませんでした……。
stable_sortにexecution::par_unseqを指定するとコンパイルできなくて悩んでましたが、今朝アップデートをかけたら普通に動きました。先にアップデート確認しとけばよかったです。
参考資料：
c++17の並列アルゴリズムライブラリ関連：
faithandbrave.hateblo.jp

tupleの独自ソート

参考資料：
minus9d.hatenablog.com

感想

c++17の新機能ってことで畳み込み式とか使ってみたかったんですが、使う場面が思いつかず。c++17の新機能は強力に見えますが、実際に使う場面が思い浮かばない辺りが何とも……研究室の他の方のプログラムに導入して改善できそうならやってみたいですね。
並列処理に関してはthreadだのOpenMPだの試してみましたが、algorithm系を除けば結局parallel_forが使いやすいかなと感じました。プラットフォームに依存しないという意味ではOpenMPの方がいいかもしれませんが……。
後Clionの構文チェックが結構いい加減というか、tupleへのgetなんかにエラー表示が出て色々とよろしくない感じがあったのがダルかったです。
そもそもc++の書き方を大半忘れていて全然書けなかったのが辛い……。
もっと良いアルゴリズムも有るようですが、力尽きたのでここまで。