Rustで作ってみよう -- HyperLogLogと並列処理で、ウィキペディア全記事のユニーク単語数を見積もる（その１）

今年の始め、私が Rust を習いはじめのころ、手本となるプログラムがあまり見つからないことが不満でした。GitHub で探せば、Rust で書かれた実用的なライブラリーが数多く見つかりますが、それらを読むのは入門者にとっては敷居が高過ぎます。私が欲しかったのは、学習用に書かれたプログラムで、入門者が手軽に試せて、いろいろといじれるプログラム例でした。

そんなわけで、そういうプログラム例を書いてみようと思います。２回に分けて、Rust で簡単なツールを作ります。

今回は乱択アルゴリズムの一種である、probability cardinarity estimatior（確率的カーディナリティ推定機）を実装します。HyperLogLog という名前のデータ構造です。

次回は HyperLogLog を使って、ウィキペディア日本語版の全記事で使われているユニーク単語数を、プラスマイナス約 2% の誤差で推定します。ちなみに、ウィキペディア日本語版には約100万個の記事があり、そのアーカイブファイルのサイズは、bzip2 形式で圧縮された状態で 2.2GB くらいあります。これを次回、マルチコアによる並列処理で、DB などは使わずにメモリー上でカウントします。また、日本語の単語の切り出しと品詞の判定には、MeCab の C API を使います。

（本文の後半がまだ書きかけですが、Rust Advent Calendar １日目の記事ですので、ここで公開して、徐々に書き進めていくことにします。最近、こんなのばっかりで、しかも、いつも未完ですみません）


# HyperLogLog とは？

HyperLogLog は確率的カーディナリティ推定機の一種です。そのアルゴリズムは、こちらの論文で解説されています。

- [Flajolet, P.; Fusy, E.; Gandouet, O.; Meunier, F. (2007) -- HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm](http://algo.inria.fr/flajolet/Publications/FlFuGaMe07.pdf)（PDF形式）

では「カーディナリティ」とはなんでしょうか？ IT、特にデータベースの世界では、この用語について微妙に異なる解釈があるようですが、cardinarity estimatior の世界に限って言うと、「集合に含まれる値の、異なり数」と考えてよさそうです。平たく言うと「distinct な値の数」です。

HyperLogLog は、確率的な手法に基づいて、与えられた値の集合のカーディナリティを推測するデータ構造です。例えば、HyperLogLog に、ウェブサイトのアクセスログのクライアント IP アドレスを与えると、ユニークビジター数がある程度の誤差の範囲で推定できます。

ところで、もし HyperLogLog のようなものがなかったとしたら、どうやって異なり数を求めるでしょうか。素朴な方法として、HashSet を使うことが考えられます。HashSet に IP アドレスのような値をどんどん詰めていくと重複が取り除かれた値が残りますので、その要素数を調べればいいわけです。この方法はシンプルですが、異なり数の分だけメモリーを消費するという欠点があります。データ量が膨大な場合、メモリーを使い果たしてしまうかもしれません。

一方、HyperLogLog は、わずか約 1.5KB のメモリー使用量で、十億を上回る数の異なり数を、プラスマイナス２％程度の誤差で推定できます。

また、HyperLogLog アルゴリズムでは、２つの HyperLogLog データ構造（以降「推定機」と呼びます）を１つにマージできます。例えば、先ほどのユニークビジターの例で、仮にアクセスログが大量にあって１台のマシンで処理しきれなかったとしましょう。こんな時は、複数台のマシンにアクセスログを分配し、それぞれに個別の推定機を持たせて処理します。処理が終わったら、全マシンから推定機を集めてマージすることで、全体のユニークビジター数が求められます。

他にも、日毎に推定機を作成しておけば、後で、月間や年間のユニークビジター数が知りたくなったときでも、マージにより瞬時に推定できます。

# 今回作るもの

では、今回どんなものができるか、その実行結果をお見せしましょう。ここでは、HyperLogLog に50万件のユニーク値をインサートしています。

```
% cargo test --release -- --nocapture --test-threads 1
...

running 2 tests
test tests::intermediate_range_with_merge_operation ...

=== Loading 500000 items to hll1 and hll2.  <= ２つの推定機それぞれに、ユニーク値を
                                                50万件ずつインサート。

（１つ目の推定機の情報）

======================================================
error rate:
  typical: 1.6250%   <== 典型的な誤差
  actual:  0.2573%   <== 実際の誤差（実行する毎に変わる）

HyperLogLog
  estimated cardinality: 498713.696162264  <== 推定値
  estimation method:     HyperLogLog       <== 推定アルゴリズム
  -----------------------------------------------------
  b:      12 bits (typical error rate: 1.625%)
  m:      4096 registers
  alpha:  0.7211100396160289
  hasher: (1527868593550391469, 4150425632492557278)

  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: **** 86
  6: *********************** 538
  7: **************************************** 963
  8: **************************************** 971
  9: **************************** 667
 10: ***************** 401
 11: ********** 238
 12: ***** 121
 13: *** 49
 14: ** 32
 15: * 18
 16: * 5
 17: * 5
 18: 0
 19: * 2

...
```

最後のヒストグラム（横棒グラフ）は、レジスターと呼ばれる「カウンターの集合体」に格納されたカウント値の分布を表しています。上の例ですと、カウント値として 0 が入っているレジスタは 0 個。1 が入っているレジスタも 0 個で、5 が入っているレジスタは 86 個、6 が入っているレジスタは 538 個でした。ここではレジスタが合計で 4096 個あり、それらに入っていた値の最大値は 19 でした。

このカウント値の意味については、次の章で説明します。


# HyperLogLog のしくみ

HyperLogLog アルゴリズムの核は以下の２つです。

- ハッシュ関数によって生成されるランダム性を持つ値
- 確率に基づく推定


## ハッシュ関数によって生成されるランダム性を持つ値

ハッシュ関数がどういうものかは、知っている人が多いと思います。あるデータ（数値でも文字列でも）を与えると、それに対応する数値（ハッシュ値）を返す関数です。

ハッシュ関数にはさまざまながアルゴリズムがあります。例えば、Rust の標準ライブラリでデフォルトとして使用されているハッシュ関数は [SipHash](https://en.wikipedia.org/wiki/SipHash) です。これは任意のバイト数のデータと、128ビットの秘密鍵を入力にして、ハッシュ値として64ビットの符号なし整数を返します。

ハッシュ関数には、そのアルゴリズムを評価するための指針がいくつかあり、代表的なものとして以下のものがあります。（ウィキペディアの「[ハッシュ関数](https://ja.wikipedia.org/wiki/%E3%83%8F%E3%83%83%E3%82%B7%E3%83%A5%E9%96%A2%E6%95%B0#.E7.89.B9.E6.80.A7)」より引用）

- **一様性**：良いハッシュ関数は、考えられる入力範囲が出力範囲全体になるべく一様に分布するようにマッピングを行う。つまり、出力範囲のそれぞれのハッシュ値はほぼ同じ確率で生成されるべきである。このような条件があるのは、異なる入力が同じハッシュ値にマッピングされてしまう「衝突」が発生すると、ハッシュに基づく各種技法のコストは衝突発生回数と共に増大するためである。

- **連続性**：（等しいデータではなく）類似するデータを探索する用途では、ハッシュ関数は可能な限り連続となっているべきである。少しだけ異なる入力に対しては、同じハッシュ値かごく近いハッシュ値を生成すべきである。なお、連続性はチェックサムや暗号学的ハッシュ関数などにとっては不適切な特性である。ハッシュ関数に連続性が必要となる用途は、線型探索を使うハッシュテーブルなどの用途である。

HyperLogLog ではハッシュ関数を用いますが、そこでは、一様性があって、連続性が **ない** ものが求められます。SipHash はその要件を満たしています。このようなハッシュ関数は、ある意味、乱数生成機に近いものとなります。類似するデータを与えた時に、非連続で、64ビット符号なし整数がとりうる値の範囲に一様に分布したハッシュ値が得られます。つまり、生成されたハッシュ値が、まるで乱数のように見えます。


## 確率に基づく推定

先ほどヒストグラムを見せました。これは、レジスターと呼ばれる「カウンターの集合体」に格納されたカウント値の分布を表しています。

```
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: **** 86
  6: *********************** 538
  7: **************************************** 963
  8: **************************************** 971
  9: **************************** 667
 10: ***************** 401
 11: ********** 238
 12: ***** 121
 13: *** 49
 14: ** 32
 15: * 18
 16: * 5
 17: * 5
 18: 0
 19: * 2
```

このカウント値は、何を意味しているのでしょうか？ 例えるなら、コインを連続して投げて、表が出るまでに投げた回数です。例えば、5 が入っているレジスタでは、コインを投げたところ、４回連続して裏が出て、５回目にやっと表が出たことを表しています。

SipHash は、ハッシュ値として64ビットの符号なし整数を返し、その値は、まるで乱数のように、ハッシュ値がとりうる値の範囲に一様に分布していると言いました。このハッシュ値を２進数で表してみます。例えば、`"test01"` という文字列と、`"test02"` という文字列について、ハッシュ値を求めたところ、以下のような２進数が得られたとします。

**ハッシュ値の例**

- `"test01"` → 0100 1001 1101 0010 ... 1011
- `"test02"` → 1001 1011 1110 0110 ... 0011

この数列を左から見て、0 はコインの裏、1 はコインの表だとします。

**コイン投げの結果**

- `"test01"` → 01 → ２回目で表が出た
- `"test02"` → 1  → １回目で表が出た

HyperLogLog では、値をインサートする度に、その値に対するハッシュ値を求め、コイン投げのように、左から何ビット目で 1 が出現したかを調べます。そして、レジスタ毎に、その最大値を記録します。

そしてここが重要なのですが、ハッシュ関数は同じ入力データに対して、同じハッシュ値を出力しますので、同じ入力データが何度現れても、その最大値の記録には影響しないことになります。これによりデータの重複がなかったものとして取り除かれ、**結果的に、異なり値の数の分だけ、コイン投げをしたことと同じになります**。

上のヒストグラムでは 19回と記録されていたレジスタもありますが、これはつまり、裏が18回も続いて出たということです。その確率が相当低いことは感覚で分かるのではないでしょうか？ 一体何回コイン投げをしたのでしょうか。裏が18回続けて出る確率は、2の18乗分の1、つまり、26万2千144分の1です。ということは、確率的に考えると、平均的に、異なり値を26万2千144種類処理すると、この記録が出ると予想できます。言い換えると、異なり値が26万2千144種類くらいあったのではないかと推定できます。

これが、HyperLogLog の基本的なしくみです。ハッシュ関数によりデータの重複を取り除きつつ乱数を得て、ビット列に1が出現する位置の最大値を記録することで、確率的な推定を行います。

ここで１つ問題があります。裏が18回続けて出る確率は26万2千144分の1だといっても、実際に何回目のハッシュ値でこれが出るかは、やってみないとわかりません。運が良ければ（運が悪ければ？）、ハッシュ値を10回求めただけでも、出てしまうかもしれません。異なり値が10個しかないのに、それを26万2千144個だと予想するのは大問題です。

これを解決するために、レジスターというカウンタの集合体が使われています。ハッシュ値の右側のいくつかのビットを元に、ハッシュ値もグループ分けします。仮に4ビットでグループ分けするとしたら、

**ハッシュ値のグループ分け**

- `"test01"` → 0100 1001 1101 0010 ... 1011 → グループ 1011（10進数なら 11）
- `"test02"` → 1001 1011 1110 0110 ... 0011 → グループ 0011（10進数なら 3）

4ビットですと16のグループに分かれますが、このグループごとに 1 が出現した位置の最大値を記録するようにして確率を求め、それらの平均を取るようにすれば、仮にどこかのグループですごく運のいいことが起こっても、推定が大きく外れることを避けられます。グループの数が増えれば増えるほど、予測される精度は上がります。

なお、この平均には、一般的に用いられている「算術平均」ではなく、「調和平均」が用いられます。調和平均は、確率の平均を求めるのに適しています。

HyperLogLog では、このグループ分けに使うビット数を b という名前で表し、その値は4から16までの範囲で指定できます。例えば、先ほどの実行例では、

```
HyperLogLog
  estimated cardinality: 498713.696162264  <== 推定値
  estimation method:     HyperLogLog       <== 推定アルゴリズム
  -----------------------------------------------------
  b:      12 bits (typical error rate: 1.625%)
  m:      4096 registers
  alpha:  0.7211100396160289
  hasher: (1527868593550391469, 4150425632492557278)
```

b として12ビットを指定し、その典型的な誤差は 1.625% となっていました。


# HyperLogLog の基本機能の実装

## プロジェクトの作成

それでは、プログラムを書いていきましょう。まず、プロジェクトを作成します。

```
$ cd 適当なディレクトリ
$ cargo new hyperloglog
$ cd hyperloglog
```

## `HyperLogLog` 構造体と `new()` メソッドを実装する

`src/lib.rs` をエディターで開いて、`HyperLogLog` 構造体を定義します。フィールドとして、論文の Fig. 3 に登場する `b`、`m`、`alpha`、`regesters` に加えて、SipHaher の初期化に使うキーなども用意します。

```rust:src/lib.rs
/// `HyperLogLog` オブジェクト
pub struct HyperLogLog {
    // レジスタのアドレッシングに使う２進数のビット数。
    // 範囲は 4以上、16以下で、大きいほど見積もり誤差が小さくなるが、その分メモリを使用する。
    b: u8,
    // usize 型のハッシュ値の右からbビットを取り出すためのマスク
    b_mask: usize,
    // レジスタの数（２のb乗）。例：b = 4 → 16、b = 16 → 65536
    m: usize,
    alpha: f64,
    // レジスタ。サイズが m バイトのバイト配列
    registers: Vec<u8>,
    // SipHasher の初期化に使うキー
    hasher_key0: u64,
    hasher_key1: u64,
}
```

Hash 関数として `SipHasher` を使用しますが、これを初期化するために乱数生成機が必要です。rand クレートを使用しましょう。Cargo.toml の `[dependencies]` セクションに rand クレートを追加します。

```toml:Cargo.toml
[dependencies]
rand = "0.3.14"
```

src/lib.rs で use しましょう。

```rust:src/lib.rs
extern crate rand;

use rand::Rng;
```

次に、この構造体を初期化するための `new()` メソッドを定義します。

論文の Fig. 3 によると、HyperLogLog には、初期化パラーメーター `b` があり、これは、レジスタのアドレッシングに使用するビット数を定めます。また、その範囲は 4 から 16 までです。まず、そこまで実装します。

```rust:src/lib.rs
use std::error::Error;

impl HyperLogLog {

    /// `HyperLogLog` オブジェクトを作成する。b で指定したビット数をレジスタの
    /// アドレッシングに使用する。b の範囲は 4以上、16以下でなければならない。
    /// 範囲外なら `Err` を返す。
    pub fn new(b: u8) -> Result<Self, Box<Error>> {
        if b < 4 || b > 16 {
            return Err(From::from(format!("b must be between 4 and 16. b = {}", b)));
        }
```

このように、b の値が範囲外の時は、`Result::Err<Box<Error>>` を返すようにしました。

次に、Fig. 3 によると、構造体のフィールド `m` は 2のb乗。シフト演算で実装しました。

```rust:src/lib.rs
    pub fn new(b: u8) -> Result<Self, Box<Error>> {
        if b < 4 || b > 16 {
            return Err(From::from(format!("b must be between 4 and 16. b = {}", b)));
        }

        let m = 1 << b;
```

次に、構造体のフィールド `alpha` は Fig.3 に書かれている定数と式使います。`get_alpha()` という別の関数にしました。

```rust:src/lib.rs
/// ビット数 b に対応する α 値を返す。
fn get_alpha(b: u8) -> Result<f64, Box<Error>> {
    if b < 4 || b > 16 {
        Err(From::from(format!("b must be between 4 and 16. b = {}", b)))
    } else {
        Ok(match b {
            4 => 0.673, // α16
            5 => 0.697, // α32
            6 => 0.709, // α64
            _ => 0.7213 / (1.0 + 1.079 / (1 << b) as f64),
        })
    }
}
```

呼び出し側は `get_alpha(b)?` というように、`?` 付きで呼び出します。これは Rust 1.13 から導入された新しい演算子で、Rust 1.0 から備わっている `try!` マクロと同じような動作をします。`get_alpha()` から返された値が `Ok<f64>` 型なら、それを unwrap して f64 値を取り出します。また、もし `get_alpha()` から返された値が `Err<Box<Error>>` 型なら、`new()` 関数から早期リターンして、呼び出し元にエラー値を返します。

```rust:src/lib.rs
impl HyperLogLog {

    pub fn new(b: u8) -> Result<Self, Box<Error>> {
        if b < 4 || b > 16 {
            return Err(From::from(format!("b must be between 4 and 16. b = {}", b)));
        }

        let m = 1 << b;
        let alpha = get_alpha(b)?;
```

最後に Hasher 初期化用の乱数を生成して、`HyperLogLog` 構造体を返します。

```rust:src/lib.rs
use std::error::Error;

impl HyperLogLog {

    /// `HyperLogLog` オブジェクトを作成する。b で指定したビット数をレジスタの
    /// アドレッシングに使用する。b の範囲は 4以上、16以下でなければならない。
    /// 範囲外なら `Err` を返す。
    pub fn new(b: u8) -> Result<Self, Box<Error>> {
        if b < 4 || b > 16 {
            return Err(From::from(format!("b must be between 4 and 16. b = {}", b)));
        }

        let m = 1 << b;
        let alpha = get_alpha(b)?;
        // hasher_key0, key1 を初期化するための乱数ジェネレータ。
        let mut rng = rand::OsRng::new().map_err(|e| format!("Failed to create an OS RNG: {}", e))?;

        Ok(HyperLogLog {
            alpha: alpha,
            b: b,
            b_mask: m - 1,
            m: m,
            registers: vec![0; m],
            hasher_key0: rng.gen(),
            hasher_key1: rng.gen(),
        })
    }
}
```

簡単なテストケースを書いて、実行してみましょう。

```rust:src/lib.rs
// テストケース
#[cfg(test)]
mod tests {

    use super::*;

    #[test]
    fn create_hll() {
        assert!(HyperLogLog::new(3).is_err());
        assert!(HyperLogLog::new(17).is_err());

        let hll = HyperLogLog::new(4);
        assert!(hll.is_ok());

        let hll = hll.unwrap();
        assert_eq!(hll.b, 4);
        assert_eq!(hll.m, 2_usize.pow(4));
        assert_eq!(hll.alpha, 0.673);
        assert_eq!(hll.registers.len(), 2_usize.pow(4));

        assert!(HyperLogLog::new(16).is_ok());
    }
}
```

```
$ cargo test
... （後ほど記入します）

```


## `insert()` メソッドを実装する

次は HyperLogLog に新しい要素を追加する `insert()` メソッドを実装しましょう。まず `SipHasher` を use します。

```rust:src/lib.rs
use std::hash::{Hash, Hasher};

// SipHasher は Rust 1.13.0 で非推奨となった。しかし、それを置き換える SipHasher24 は
// 現状では非安定(unstable) なため、安定版の Rust リリースでは使用できない。
#[allow(deprecated)]
use std::hash::SipHasher;
```

論文の Fig. 3 に沿って実装していきます。値 `value` の hash 値を変数 `x` に束縛。`b_mask` フィールドの値を使って、`x` の値の右から b ビットを切り出し、これをレジスタ Vec の要素番号として使われる変数 `j` に束縛します。また、`x` の残りのビット（64ビット - b ビット）を変数 `w` に束縛します。

```rust:src/lib.rs
impl HyperLogLog {

    // ...

    /// 要素を追加する。要素は `std::hash::Hash` トレイトを実装していなければならない。
    pub fn insert<H: Hash>(&mut self, value: &H) {
        let x = self.hash(value);
        let j = x as usize & self.b_mask;
        let w = x >> self.b;


    }

    /// 与えられた value に対する 64 ビットのハッシュ値を求める。
    #[allow(deprecated)]  // SipHasher が Rust 1.13 で非推奨(deprecated) のため。
    fn hash<H: Hash>(&self, value: &H) -> u64 {
        let mut hasher = SipHasher::new_with_keys(self.hasher_key0, self.hasher_key1);
        value.hash(&mut hasher);
        hasher.finish()
    }

}
```

次に、`w` の値について、左側から何桁目に `1` が立っているビットが現れたのかを求めます。これは、別関数に切り出しましょう。

```rust:src/lib.rs
impl HyperLogLog {

    // ...

    /// 要素を追加する。要素は `std::hash::Hash` トレイトを実装していなければならない。
    pub fn insert<H: Hash>(&mut self, value: &H) {
        let x = self.hash(value);
        let j = x as usize & self.b_mask;
        let w = x >> self.b;

        let p1 = position_of_leftmost_one_bit(w, 64 - self.b);


    }

}
```

以下のように実装しました。

```rust:src/lib.rs
/// ハッシュ値（64ビット符号なし２進数）の左端から見て最初に出現した 1 の位置を返す。
/// 例：10000... → 1、00010... → 4
fn position_of_leftmost_one_bit(s: u64, max_width: u8) -> u8 {
    count_leading_zeros(s, max_width) + 1
}

/// ハッシュ値（64ビット符号なし２進数）左端に連続して並んでいる 0 の個数を返す。
/// 例：10000... → 0、00010... → 3
fn count_leading_zeros(mut s: u64, max_width: u8) -> u8 {
    let mut lz = max_width;
    while s != 0 {
        lz -= 1;
        s >>= 1;
    }
    lz
}
```

最後に、`position_of_leftmost_one_bit()` が返した値と、`registers[j]` の値を比較して、大きい方の値を `registers[j]` に格納します。

```rust:src/lib.rs
impl HyperLogLog {

    // ...

    /// 要素を追加する。要素は `std::hash::Hash` トレイトを実装していなければならない。
    pub fn insert<H: Hash>(&mut self, value: &H) {
        let x = self.hash(value);
        let j = x as usize & self.b_mask;
        let w = x >> self.b;

        let p1 = position_of_leftmost_one_bit(w, 64 - self.b);
        let p2 = &mut self.registers[j];
        if *p2 < p1 {
            *p2 = p1;
        }
    }

}
```

テストケースを追加しましょう。

```rust:src/lib.rs
mod tests {

    use super::*;

    #[test]
    fn small_range() {
        let mut hll = HyperLogLog::new(12).unwrap();
        let items = ["test1", "test2", "test3", "test2", "test2", "test2"];

        println!("\n=== Loading {} items.\n", items.len());
        for item in &items {
            hll.insert(item);
        }
    }
}
```

## `cardinality()` メソッドを実装する

カーディナリティを推定する `cardinality()` メソッドを実装しましょう。Fig. 3 に沿って実装します。

```rust:src/lib.rs
/// 推定アルゴリズム。デバッグ出力用
#[derive(Debug)]
pub enum Estimator {
    HyperLogLog,
    LinearCounting, // スモールレンジの見積もりに使用する。
}

impl HyperLogLog {

    // ..

    /// カーディナリティの見積もり値を返す。
    pub fn cardinality(&self) -> f64 {
        estimate_cardinality(self).0
    }

}

// ...

/// カーディナリティを推定し、その値と、見積もりに使用したアルゴリズムを返す。
/// スモールレンジでは `Linear Counting` アルゴリズムを使用し、それを超えるレンジでは
/// `HyperLogLog` アルゴリズムを使用する。ここまでは論文の通り。
/// しかし、論文にあるラージレンジ補正は行なわない。なぜなら、本実装では、32 ビットの
/// ハッシュ値の代わりに 64 ビットのハッシュ値を使用しており、ハッシュ値が衝突する頻度が
/// 極めて低いと予想されるため。
fn estimate_cardinality(hll: &HyperLogLog) -> (f64, Estimator) {
    let m_f64 = hll.m as f64;
    // まず `HyperLogLog` アルゴリズムによる見積もり値を算出する。
    let est = raw_hyperloglog_estimate(hll.alpha, m_f64, &hll.registers);

    if est < (5.0 / 2.0 * m_f64) {
        // スモールレンジの見積もりを行う。もし値が 0 のレジスタが１つでもあるなら、
        // `Linear Counting` アルゴリズムで見積もり直す。
        match count_zero_registers(&hll.registers) {
            0 => (est, Estimator::HyperLogLog),
            v => (linear_counting_estimate(m_f64, v as f64), Estimator::LinearCounting),
        }
    } else {
        (est, Estimator::HyperLogLog)
    }
}

/// 値が 0 のレジスタの個数を返す。
fn count_zero_registers(registers: &[u8]) -> usize {
    registers.iter().filter(|&x| *x == 0).count()
}

```

個々の推定アルゴリズムの実装は以下の通り。

```rust:src/lib.rs
/// `HyperLogLog` アルゴリズムによる未補正の見積もり値を算出する。
fn raw_hyperloglog_estimate(alpha: f64, m: f64, registers: &[u8]) -> f64 {
    let sum = registers.iter().map(|&x| 2.0f64.powi(-(x as i32))).sum::<f64>();
    alpha * m * m / sum
}

/// `Linear Counting` アルゴリズムによる見積もり値を算出する。
fn linear_counting_estimate(m: f64, number_of_zero_registers: f64) -> f64 {
    m * (m / number_of_zero_registers).ln()
}
```

テストケースを追加します（後ほど追加します）


# デバッグ出力やヒストグラムを追加する

## デバッグ出力を追加する

```rust:src/lib.rs
    /// b から予想される典型的なエラー率を返す。
    pub fn typical_error_rate(&self) -> f64 {
        1.04 / (self.m as f64).sqrt()
    }
```

```rust:src/lib.rs
use std::fmt;
```

```rust:src/lib.rs
/// `HyperLogLog` のデバッグ出力用文字列を返す。
impl fmt::Debug for HyperLogLog {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let (est, est_method) = estimate_cardinality(self);
        write!(f,
               r#"HyperLogLog
  estimated cardinality: {}
  estimation method:     {:?}
  -----------------------------------------------------
  b:      {} bits (typical error rate: {}%)
  m:      {} registers
  alpha:  {}
  hasher: ({}, {})"#,
               est,
               est_method,
               self.b,
               self.typical_error_rate() * 100.0,
               self.m,
               self.alpha,
               self.hasher_key0,
               self.hasher_key1)
    }
}
```


## ヒストグラムを描けるようにする

```rust:src/lib.rs
use std::collections::BTreeMap;
```

```rust:src/lib.rs
impl HyperLogLog {

    // ...

    /// レジスタに格納された値について、その分布を示すヒストグラムを返す。
    pub fn histgram_of_register_value_distribution(&self) -> String {
        let mut histgram = Vec::new();

        let mut map = BTreeMap::new();
        for x in &self.registers {
            let count = map.entry(*x).or_insert(0);
            *count += 1;
        }

        if let (Some(last_reg_value), Some(max_count)) = (map.keys().last(), map.values().max()) {
            // グラフの最大幅 = 40 文字
            let width = 40.0;
            let rate = width / (*max_count as f64);

            for i in 0..(last_reg_value + 1) {
                let mut line = format!("{:3}: ", i);

                if let Some(count) = map.get(&i) {
                    // アスタリスク(*) で横棒を描く
                    let h_bar = std::iter::repeat("*")
                        .take((*count as f64 * rate).ceil() as usize)
                        .collect::<String>();
                    line.push_str(&h_bar);
                    line.push_str(&format!(" {}", count));
                } else {
                    line.push_str("0");
                };

                histgram.push(line);
            }
        }
        histgram.join("\n")
    }
}
```

# `merge()` メソッドを実装する

```rust:src/lib.rs
impl HyperLogLog {

    // ...

    /// `template` で与えられた既存の `HyperLogLog` オブジェクトと同じ設定を持つ
    /// `HyperLogLog` オブジェクトを作成する。
    pub fn from_template(template: &HyperLogLog) -> Self {
        let m = template.m;
        HyperLogLog {
            alpha: template.alpha,
            b: template.b,
            b_mask: m - 1,
            m: m,
            registers: vec![0; m],
            hasher_key0: template.hasher_key0,
            hasher_key1: template.hasher_key1,
        }
    }

    /// `self` で示される `HyperLogLog` オブジェクトへ、`other` で示される別の
    /// `HyperLogLog` オブジェクトをマージする。両オブジェクトの設定が異なる場合は
    /// `Err` を返す。
    pub fn merge(&mut self, other: &HyperLogLog) -> Result<(), Box<Error>> {
        if self.b == other.b && self.m == other.m && self.hasher_key0 == other.hasher_key0 &&
           self.hasher_key1 == other.hasher_key1 {
            for (p1, p2) in self.registers.iter_mut().zip(other.registers.iter()) {
                if *p1 < *p2 {
                    *p1 = *p2
                }
            }
            Ok(())
        } else {
            Err(From::from(format!("Specs does not match. \
                                    b: {}|{}, m: {}|{}, hasher: ({},{})|({},{})",
                                   self.b,
                                   other.b,
                                   self.m,
                                   other.m,
                                   self.hasher_key0,
                                   self.hasher_key1,
                                   other.hasher_key0,
                                   other.hasher_key1)))
        }
    }

}
```

# まとめ

- HyperLogLog の基本的な仕組みを説明し、Rust で実装してみました。
- 単にカーディナリティを推測できるだけでなく、ヒストグラムなどデバッグ用の出力を充実させ、内部の動作が学習できるようにしました。
- 次回は HyperLogLog を使って、ウィキペディア日本語版の全記事で使われているユニーク単語数を、プラスマイナス約 2% の誤差で推定します。