C++ による数値シミュレーション入門 (4/4)｜実装例

東工大情報理工学院 高安研究室 で開催されている新入生向けのプログラミングゼミの資料を一部公開します．

本稿は C++ による数値シミュレーション入門 の第 4 部です．

• (1/4) 導入
• (2/4) C++ の基本構文
• (3/4) 標準ライブラリ
• (4/4) 実装例 ← 今回

ここでは，次の2つのモデルの実装を行ないます．

1. Barabashi-Albert モデル
2. ディーラーモデル（課題）

Barabashi-Albert モデル

現実のネットワークの生成モデルである Barabashi-Albert モデル（以下，BA モデル）を簡単に紹介した後，その実装例を step by step で紹介します．

[5] Barabási, Albert-László, and Réka Albert. "Emergence of scaling in random networks." science 286.5439 (1999): 509-512. https://www.science.org/doi/10.1126/science.286.5439.509

現実のネットワークのスケールフリー性

ネットワークはノードとノードがリンクによって結ばれてできています．例えば，正方格子は全てのノードが一様に，同じ数のノードとリンクで繋がれています．しかし，現実のネットワーク（例えば，SNS のフォロー関係のネットワークやタンパク質の相互作用のネットワーク）にはそのような均質性はなく，1つのノードが大部分のリンクを占有しています．このような不均質なつながりはスケールフリー性と呼ばれ，ネットワークの複雑さを高めています．スケールフリー性をより正確に表すと，それぞれのノードのリンクでつながったノードの個数を $k$ としたとき，その分布が

$$P(>k) \sim k^{-\alpha}$$

のように，べき分布で記述される性質のことです．$k$ はノードの次数 (degree) とも呼ばれます．次数はノードの規模を表す一つの指標になります．正の実数 $\alpha$ はべき指数と呼ばれ，$\alpha$ が小さいほど次数の偏りが大きくなります．スケールフリー性をもつ現実のネットワークは，$1<\alpha<3$ の値を取ることが多いです．

BA モデル

BA モデルは，スケールフリー性をもつネットワークが生成される過程を記述した，1番最初のモデルです．スケールフリー性は，富めるものがさらに豊かになる性質である，Preferential Attachment (優先的接続) をネットワークの成長過程に導入しました．Twitter を例に挙げれば，フォロワー数の多い人ほど（有名人ほど）フォローされやすいという性質です．具体的には次の通りです．

1. $N_0$ 個のノードからなる（小さな）ネットワークを用意する．すべてのノードがリンクを介してつながっていれば形状は何でも良い．リングでも，すべてのノード同士がつながった完全グラフでも．
2. $N$ 個のノードを順番にネットワークに追加する．
   1. ノードは $m$ 本のリンクをもつ（この時点ではリンクの先はどのノードにもつながっていない）
   2. それぞれのリンクの接続先を，ネットワーク上にあるノードから次の確率で決める．
      $$\Pi_i = \frac{k_i}{\sum_j k_j} \propto k_i$$
      すなわち，次数の大きさに比例する確率でリンクをつなぐ先を選ぶ．このルールは，次数が大きいノードを優先的に選ぶことから，Preferential Attachment (優先的接続) と呼ばれる．
3. 十分大きな $N$ を取ると，$N_0$ によらず，$P(>k) \sim k^{-2}$ のスケールフリーネットワークが生成される．

このモデルによって生成されるネットワークが $\alpha=2$ の次数分布をもつことは，ネットワーク科学の参考書などで勉強してください．また，べき指数が $\alpha = 2$ でないネットワークは，Preferential Attachment の確率の与え方を

$$\Pi_i \propto k_i + c$$

$$\Pi_i \propto k_i^c$$

などと変えることで実現できます．

それでは，BA モデルを実装していきます．考えるべき点は以下です．

1. ネットワークのデータ構造
2. $N_0$ 個のノードからなる初期ネットワークの生成
3. Preferential Attachment の実装
4. （任意）ノード追加時に，自己ループができないようにリンクの接続先を選ぶ
5. （任意）ノード追加時に，多重リンクができないようにリンクの接続先を選ぶ
6. 生成したネットワークをファイルに出力する

ネットワークのデータ構造

ネットワークの構成要素であるノードは，ノード数が $N$ であれば，0から$N-1$ の連番で識別することにします．ネットワークはノード同士を結ぶリンクの集合ですから，データ構造は

// 例: 0-1-2
links = [(0,1),(1,2),(2,3)]

のような組の配列で記述できそうです．リンクに向きがない場合は，一方向のみを含める場合と，両方向を一回ずつ含める場合がありますが，ここでは前者を採用することにします．このとき，数字の小さいノードを先に書くなどのルールを決めておくと良いでしょう．

また，偶数番目とその次を1組のリンクだと思うと，() を取り払って

// 例: 0-1-2
links = [0,1,1,2,2,3]

のようにデータ構造を定義できそうです．$L$ 本のリンクが $2L$ 個の要素からなる配列で表現できるわけです．今回はこのデータ構造を採用して，ネットワークを次で定義することにします．

// 例: 0-1-2
std::vector<int> links = {0,1,1,2,2,3};

リンクを追加していく必要があるので，可変長配列 std::vector を使っています．

初期ネットワークの生成

$N_0$ 個の初期ネットワークとして，完全リンクを作ることにします．例えば，$N_0=4$ の場合は

(0,1),(0,2),(0,3),
      (1,2),(1,3),
            (2,3)

の合計6本のリンクが作られます．これは二重 for 文を使って簡単に書けそうです．

// Assume undirected graph
// Node labels: 0,...,N-1
std::vector<int> generate_complete_graph(int N)
{
    std::vector<int> links;
    for (int i = 0; i < N; ++i)
    {
        for (int j = i + 1; j < N; ++j)
        {
            links.push_back(i);
            links.push_back(j);
        }
    }
    return links;
}

links に (i,j) リンクを追加することは，この後も何度も出てきそうですので，次のように関数化してしまっても良いでしょう．

void add_link(std::vector<int> &links, int i, int j)
{
    // add i -> j link (i < j)
    links.push_back(i);
    links.push_back(j);
}

関数には links を参照渡しする必要があります．& を付けず値渡ししてしまうと，コピーされた links に対して(i,j) リンクが追加されてしまいます．

Preferential Attachement

ノードを追加する際，リンクの先をネットワーク上のノードから次数に比例する確率

$$\Pi_i = \frac{k_i}{\sum_j k_j} \propto k_i$$

で選ぶことをそのままプログラムで再現しようとすると，次数の配列を用意して，離散分布を利用する

std::vector<int> degrees = {k0, k1, ...};
std::discrete_distribution<int> discrete(degrees);
int target = discrete(gen);  // gen は乱数生成器

といった発想になるかと思います．ノードを追加する度に degrees をするのが面倒ですが，いちおう実装はできそうです．ただ，もっと簡単な方法があります．

実は，links から一様ランダムに要素を取り出すだけで preferential attachment を実装できます．links にはノード $i$ がちょうど $k_i$ 回出てきます．次数の和は $\sum_i k_i = 2L$ であり，ちょうど links.size() に一致します．つまり，links からランダムに要素を選んだ時に，ノード $i$ が選ばれる確率は

$$\frac{k_i}{2L} = \frac{k_i}{\sum_j k_j} = \Pi_i$$

に一致します．このことを利用すると，preferential attachment による接続先の選択は次のように書けます．

std::uniform_int_distribution<int> uni(0, links.size()-1);
int target = links[uni(gen)];  // gen は乱数生成器

自己ループを避ける（任意）

新規ノードの接続先の数 $m$ が 2 以上のとき，実装方法によっては自己ループが発生します．

• 「接続先を 1 つ選んでリンクを追加することを $m$ 回繰り返す」という実装だと，2 個目以降の接続先として自身が選ばれる可能性があります．
• 「$m$ 個の接続先をすべて選んでから $m$ 本のリンクを追加する」という順番にすると，接続先として自身が選ばれることはありません．

多重リンクを避ける（任意）

$m$ 個の接続先を選ぶ先，接続先に重複があると多重リンクができてしまう．$m$ 個の接続先を順に決定する際，「新しく選んだものが既に選ばれているか判定し，選ばれていたら選び直す」という機構が必要です．判定用に次のような contain 関数を用意しておくと便利です．

// Determine if a vector contains specified target
bool contains(std::vector<int> &vec, int target)
{
    for (auto element : vec)
    {
        if (element == target)
            return true;
    }
    return false;
}

この C++ のコードは，標準ライブラリの std::find を使うことでもっとシンプルに書くこともできます．

#include <vector>
#include <algorithm> // std::find を使うために必要

bool contains(const std::vector<int> &vec, int target)
{
    return std::find(vec.begin(), vec.end(), target) != vec.end();
}

ファイル出力

BA モデルによって生成されたネットワーク（リンクの集合）をファイルに出力します．ファイル出力には fstream ライブラリの std::ofstream を使います．下記のように関数化しておくと main 関数の見通しが良くなります．

void output_links(std::vector<int> &links, std::string path = "")
{
    if (path == "")
    {
        for (int i = 0; i < links.size() / 2; ++i)
            std::cout << links[2 * i] << "," << links[2 * i + 1] << std::endl;
        return;
    }
    std::ofstream ofs(path);
    for (int i = 0; i < links.size() / 2; ++i)
        ofs << links[2 * i] << "," << links[2 * i + 1] << std::endl;
}

プログラムの全体像

以上をまとめると，BA モデルは下記プログラムのように実装できます．

ba.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <random>

void add_link(std::vector<int> &links, int i, int j)
{
    // add i -> j link (i < j)
    links.push_back(i);
    links.push_back(j);
}

// Assume undirected graph
// Node labels: 0,...,N-1
std::vector<int> generate_complete_graph(int N)
{
    std::vector<int> links;
    for (int i = 0; i < N; ++i)
    {
        for (int j = i + 1; j < N; ++j)
        {
            add_link(links, i, j);
        }
    }
    return links;
}

// Determine if a vector contains specified target
bool contains(std::vector<int> &vec, int target)
{
    for (auto element : vec)
    {
        if (element == target)
            return true;
    }
    return false;
}

// BA model: generate a network (list of links)
std::vector<int> BA(std::mt19937 &gen, int N, int N0, int m)
{
    // Initialize a graph
    std::vector<int> links = generate_complete_graph(N0);

    // For N-N0 steps, add a node with m links
    // Labels of new nodoes: N0,...,N-1
    // Choose targets of m links by preferential attachment
    std::uniform_real_distribution<double> uni(0.0, 1.0);
    for (int new_node = N0; new_node < N; ++new_node)
    {
        // Select targets to connect m links 
        // by preferential attachment
        std::vector<int> targets;
        int L = links.size();
        while (targets.size() < m)
        {
            int index = int(uni(gen) * L);
            int new_target = links[index];
            // Avoid multi-links
            if (contains(targets, new_target))
                continue;
            targets.push_back(new_target);
        }
        // add m links to the network
        for (auto target: targets)
        {
            add_link(links, target, new_node);
        }
    }
    return links;
}

// output links as csv
void output_links(std::vector<int> &links, std::string path = "")
{
    if (path == "")
    {
        for (int i = 0; i < links.size() / 2; ++i)
            std::cout << links[2 * i] << "," << links[2 * i + 1] << std::endl;
        return;
    }
    std::ofstream ofs(path);
    for (int i = 0; i < links.size() / 2; ++i)
        ofs << links[2 * i] << "," << links[2 * i + 1] << std::endl;
}

int main(int argc, char *argv[])
{
    if (argc < 2)
    {
        std::cout << "Required an output directory." << std::endl;
        std::cout << "[Usage] <bin> <out_dir>" << std::endl;
        return 0;
    }
    // 乱数列の出力ファイル
    // std::string out_dir = "../result/";
    std::string out_dir = argv[1];
    std::string out_file = out_dir + "ba.csv";

    int seed = 1;
    std::mt19937 gen(seed);

    int N = 1e6;
    int N0 = 5;
    int m = 3; // 新規参入するノードのリンク数
    std::vector<int> links = BA(gen, N, N0, m);

    output_links(links, out_file);

    return 0;
}

課題 4（ディーラーモデルの実装）

ディーラーモデルを実装してください．

Yamada, Kenta, Hideki Takayasu, Takatoshi Ito, and Misako Takayasu. 2009. “Solvable Stochastic Dealer Models for Financial Markets.” Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 79 (5 Pt 1): 051120. http://dx.doi.org/10.1103/PhysRevE.79.051120

実装例：C++ によるディーラーモデルの実装