C++17で基本統計量を計算してみる

こんにちは id:mitsutaka-takada です。

今日はC++17で導入されるRanges TS(Technical Specification)ライブラリの紹介をしたいと思います。

http://en.cppreference.com/w/cpp/experimental/ranges

Rangesライブラリが導入されると、C++のプログラミングにどのような影響があるか数式をプログラミングしてみることで確認していきたいと思います。

プログラミングのお題は以下のエントリで紹介されているものを使い、Rangesライブラリ以前と以降で書き方の違いが明確になれば良いかと思います。

qiita.com

qiita.com

最大・最小値

以下が0から19までの整数から最大値(19)と最小値(0)を求めるプログラムです。

ranges::view::iotaで0から19までの整数の列をlazyに生成して、ranges::minmaxで最大値と最小値を取り出します。

#include <range/v3/view/iota.hpp>
#include <range/v3/algorithm/minmax.hpp>

#include <cassert>

int main(int argc, char *argv[])
{
    // [0,20) の整数。
    auto const intFrom0To20 = ranges::view::iota(0, 20);
    auto const minMax = ranges::minmax(intFrom0To20);

    assert(minMax.first == 0 && minMax.second == 19);
    return 0;
}

総和

ranges::view::iotaはlazyな列なので、終点を指定せずに無限長の列(iota(1))を表現することができます。

無限長の列の総和はとれないので、列の長さを制限する必要があります。LINQやHaskellのようにtakeで列の長さを指定すると有限長の列になります。

#include <range/v3/view/iota.hpp>
#include <range/v3/view/take.hpp>
#include <range/v3/numeric/accumulate.hpp>

#include <cassert>

int main(int argc, char *argv[])
{
    // 1 .. 100 までの整数。
    auto const intsFrom1To100 = ranges::view::iota(1) | ranges::view::take(100);
    // 0は合計の初期値。
    assert(ranges::accumulate(intsFrom1To100, 0)  == (1 + 100)*100/2);
    return 0;
}

平均

平均を計算する際は従来のSTLと同じようにaccumurateを使用します。

#include <range/v3/view/iota.hpp>
#include <range/v3/numeric/accumulate.hpp>

#include <cassert>

int main(int argc, char *argv[])
{
    // 1 .. 100 までの整数。
    auto const intsFrom1To100 = ranges::view::iota(1, 101);
    assert(ranges::accumulate(intsFrom1To100, 0.0)/ranges::size(intsFrom1To100)  == 50.5);
    return 0;
}

分散

平均と内積から分散を求めます。内積はranges::inner_productで計算できます。

今回の計算では整数では浮動小数点を使用しています。ranges::view::iotaは整数しかサポートしていないのでranges::view::generate_nを使用して1から100までの浮動小数点を生成してranges::to_vectorでstd::vectorに変換します。

#include <range/v3/view/generate_n.hpp>
#include <range/v3/numeric/accumulate.hpp>
#include <range/v3/numeric/inner_product.hpp>

#include <iostream>
#include <vector>

int main(int argc, char *argv[])
{
    auto const from1To100 = ranges::view::generate_n(
            [v = 1.0]() mutable { 
                auto const tmp = v;
                ++v;
                return tmp;
            },
            100) | ranges::to_vector;

    auto const avg = ranges::accumulate(from1To100, 0.0)/from1To100.size();
    auto const variance = ranges::inner_product(from1To100, from1To100, 0)/from1To100.size()- avg*avg;
    std::cout << "variance = " << variance << std::endl;
    return 0;
}

中央値

中央値はranges::n_th_elementで求めます。n_th_elementは、対象の範囲をソートした時、指定したn番目の位置にn番目の値が来るように範囲をソートするアルゴリズムです。n番目より前の要素はn番目の要素より小さいことが保証されているだけでソートはされていません。例えば、1から5までの整数がランダムに配置された範囲をn_th_elementにn=3を指定してソートするとn番目の位置には3が来ますが、１番目と２番目の要素はどちらが1でどちらが2かは保証されていません。[1, 2, 3, X, Y]かもしれませんし、[2, 1, 3, X, Y]かもしれません。STLの中でも使われる機会が少ないアルゴリズムでしょうか。

また中央値を求める前にranges::action::shuffleで範囲をランダムにしています。action下のアルゴリズムはコンテナの要素をinplaceで処理するアルゴリズムです。この場合from1To100というvectorの要素を並べかえるのでactionとなります。一方view下のアルゴリズムはコンテナ中の要素を変更しません。

#include <range/v3/view/generate_n.hpp>
#include <range/v3/iterator_range.hpp>
#include <range/v3/numeric/accumulate.hpp>
#include <range/v3/numeric/inner_product.hpp>
#include <range/v3/algorithm/nth_element.hpp>
#include <range/v3/action/shuffle.hpp>

#include <random>
#include <iostream>
#include <vector>

int main(int argc, char *argv[])
{
    auto from1To100 = ranges::view::generate_n(
            [v = 1.0]() mutable { 
                auto const tmp = v;
                ++v;
                return tmp;
            },
            100) | ranges::to_vector;
    std::mt19937 rbg;
    from1To100 = std::move(from1To100) | ranges::action::shuffle(rbg);

    double median = 0.0;
    if(from1To100.size()%2 == 0){ 
        auto const middlePoint = from1To100.begin() + from1To100.size()/2 - 1;
        ranges::nth_element(from1To100, middlePoint);
        ranges::nth_element(middlePoint + 1, middlePoint + 1, from1To100.end());
        median = (*middlePoint + *(middlePoint + 1))/2;
    } else {
        auto const middlePoint = from1To100.begin() + from1To100.size()/2;
        ranges::nth_element(from1To100, middlePoint);
        median = *middlePoint;
    }

    std::cout << " median = " << median << std::endl;
    return 0;
}

最頻値

最頻値は、まず範囲をソート(ranges::action::sort)します。その後、値が同じ者同士を１つのグループにします(ranges::view::group_by)。
同値グループをグループの長さとそのグループの値の組に変換して(ranges::view::transform)、最後にグループの長さが最大(ranges::max)のグループを選択します。

#include <range/v3/action/sort.hpp>
#include <range/v3/algorithm/max.hpp>
#include <range/v3/view/group_by.hpp>
#include <range/v3/view/transform.hpp>

#include <iostream>
#include <vector>

int main(int argc, char *argv[])
{
    struct valAndCount{
        int value;
        long count;
    };

    auto const sorted = std::vector<int>{0, 0, 1, 10, 2, 5, 1, 1} | ranges::action::sort;
    auto mode = ranges::max(
        sorted
        | ranges::view::group_by([](auto x, auto y){ return x == y;})
        | ranges::view::transform(
            [](auto const& vs){ return valAndCount{ *vs.begin(), ranges::distance(vs)}; }),
        [](auto const& lvc, auto const& rvc){
             return lvc.count < rvc.count;
         }
        );
    std::cout << mode.value << std::endl;

    return 0;
}

最後に

Range TSを利用したプログラミングはいかがでしたでしょうか。従来のSTLで見慣れた名前が多かったかと思います。

Range TSは従来のイテレータ・ベースのアルゴリズムをRangeというコンセプト（例えば、イテレータのペアもRangeです）に拡張したもの＋αなライブラリです。アルゴリズムを使う際、イテレータのペアを自分で管理しなくてよかったり、パイプ(|演算子)で複数の処理を合成していくことが可能な便利なライブラリです。

Range TSの理論的な部分は以下の本を

www.amazon.co.jp

Range TSのポテンシャルの高さについて知りたい人は、CPPCON 2015のRange TSライブラリ作者の発表がオススメです。

www.youtube.com