炭化水素を列挙してみた（多重結合を含む）

プログラミング 数学

この記事は以下の記事の続編です。もし読んでいない方は先に前回をお読みください。
smooth-pudding.hatenablog.com

前回は、単結合のみを含む炭化水素を列挙する問題に取り組んでみました。今回は、これらの炭化水素に脱水素を施して、二重結合や三重結合を含むようなものについても列挙を実行します。

今回もこちらの書籍『炭素と水素による構造式一覧』を大いに参考にしました。
techbookfest.org

例によってコードは Rust で書いており、github に上がっています。とにかくコードを見たいという方は以下からどうぞ。この記事はこのコードの解説です。
github.com

※本文中に引用したコードは最新版と異なる場合があります。

二重・三重結合の扱い

前回の記事では単結合のみを扱いました。単結合のみの構造の場合は、各要素が 0 または 1 からなる行列で表現できるというお話をしました。

例えばこの構造は

こんな行列で表されるのでした。

\begin{pmatrix}0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \end{pmatrix}

さて、先ほどとほぼ同じだけど、一か所だけ二重結合になったこんな構造を考えてみます。

これも行列で表現できるでしょうか？

答えは Yes です。単に「二重結合のエッジは2で表現する」とするだけです。具体的にはこんな行列になります。ちょうど、1 と 3 の間の結合の部分だけ 2 に置き換わった行列ですね。

\begin{pmatrix}0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 2 \\ 0 & 0 & 0 & 1 \\ 1 & 2 & 1 & 0 \end{pmatrix}

同様にして、三重結合は 3 で表現することにしましょう。四重結合はないので*1、行列の成分を 0, 1, 2, 3 に拡張するだけでよさそうです。

ということで、今回の列挙の対象は以下の条件となります。

  • 対称行列
  • 各成分は 0, 1, 2, 3 のいずれか
  • 対角成分は 0
  • 各行（および各列）の和は4以下
  • ノードの番号の対応を入れ替えて行列が一致するものは同じものとみなす
  • 連結なグラフを表現している

要素が 0/1 から 0/1/2/3 になったので、1ビットから2ビットにすればよさそうですね。コードでは matrix.rs の SymmetricTwoBitsMatrix<N>という構造体でこれを表現しています。

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct SymmetricTwoBitsMatrix<const N: usize> {
    rows: [u32; N],
}

脱水素化ストリーム

問題が広くなったのでまた頑張ってイチから全列挙・・・する必要はありません。『炭素と水素による構造式一覧』では脱水素化ストリームと呼ばれている方法により、単結合のみの全列挙の結果を活用して少ない計算で列挙することができます。

先ほど二重結合を含む分子構造の例を挙げました。この構造は単結合のみからなる構造とよく似ていて、一か所だけが二重結合になったものでした。言い方を変えれば、「単結合のみからなる構造から、二重結合も含まれる構造を作る操作をした」と言えます。

ここでひらめきます。「二重結合や三重結合が含まれている構造を列挙するのであれば、まずは全ての結合を単結合に置き換えたものを列挙して、そのあとに二重結合や三重結合に置き換えたものを作ればいいんじゃないか？」これが脱水素ストリームの肝となる考え方です。

ちなみに、先ほどの例だと操作の前の分子式が C4H8 で操作後が C4H6 となっています。水素を減らす操作ということで脱水素という名前がぴったりですね。

脱水素化ストリームを隣接行列の言葉で表現すれば、以下のようなステップになります。

  1. 単結合のみからなる構造のうちのひとつを選び、隣接行列を作る
  2. 1 になっている要素の数だけ、それを +1 した行列を作る
  3. 作った行列たちのうち、ノードの番号付けを入れ替えても一致するものは省く
  4. 残されたもののそれぞれの隣接行列に着目する
  5. 1 または 2 になっている要素の数だけ、それを +1 した行列を作る（以下繰り返し）
  6. 脱水素する対象がなくなったら終了

ここで +1 する際には以下に気を付けるようにします。

  • +1 した結果、各行の和が 4 を越えないようにする
  • (上記の条件から自動的に満たされますが*2 ) どの要素も3を超えないようにする

コードでは dehydrogenate.rs の generate_all_dehydrogenated<N> という関数で行っています*3。

pub fn generate_all_dehydrogenated<const N: usize>(
    base_matrix: SymmetricTwoBitsMatrix<N>,
    symmetry: &[[usize; N]],
) -> Vec<SymmetricTwoBitsMatrix<N>> {
    let mut all_mats = vec![base_matrix];
    let mut seen = vec![[base_matrix].iter().copied().collect::<FxHashSet<_>>()];
    let mut current = vec![base_matrix];
    let mut next = Vec::new();

    while !current.is_empty() {
        for mat in current {
            for (row, col) in mat.dehydrogenatable_bonds() {
                let mat_new = mat.create_dehydrogenated_unchecked(row, col);
                if seen
                    .iter()
                    .any(|s: &FxHashSet<SymmetricTwoBitsMatrix<N>>| s.contains(&mat_new))
                {
                    continue;
                }
                all_mats.push(mat_new);
                next.push(mat_new);

                let variants = make_variants(&mat_new, symmetry);
                seen.push(variants);
            }
        }

        current = next;
        next = Vec::new();
        seen.clear();
    }
    all_mats
}

SymmetricTwoBitsMatrix<N> の create_dehydrogenated_unchecked で実際に要素に +1 する操作を行っています。

pub fn create_dehydrogenated_unchecked(&self, row: usize, col: usize) -> Self {
    let mut new_matrix = *self;
    new_matrix.rows[row] += 1 << (2 * col);
    new_matrix.rows[col] += 1 << (2 * row);
    new_matrix
}

計算結果

こんな感じで、炭素数ごと/水素の数ごとに構造の数を表示させてみました。すべて『炭素と水素による構造式一覧』にある表と一致しました。やったぜ。

===== [C =  1] =====
#H: # of structures
 4: 1
===== [C =  2] =====
#H: # of structures
 2: 1
 4: 1
 6: 1
===== [C =  3] =====

// ~~~ (中略) ~~~

===== [C = 10] =====
#H: # of structures
 0: 4330
 2: 64352
 4: 241297
 6: 439373
 8: 488125
10: 369067
12: 201578
14: 81909
16: 24938
18: 5568
20: 852
22: 75

単結合の構造の列挙から始めていますが、計算時間はなんと6秒弱になりました*4。速い！

計算の環境は前回と同じなので省略します。前回の記事を参照してください。

高速化の工夫

せっかく Rust で書いたので、今回は自分の知識の範囲で可能な限りギリギリまで高速化してみました。その過程で行った工夫を最後にいくつか紹介します。

対称性の計算の共通化

脱水素化したものたちで異なる構造だけ残すために、ノードの番号付けを入れ替える操作を行う必要があります。N! 通り試すと大変ですが、よく考えてみると、単結合の隣接行列を変えてしまうような番号の入れ替えはそもそも考えなくてよいことが分かります。仮にそのようなものを考慮しても、脱水素化した別の候補とは絶対に一致しないからです。

そのため、脱水素化ストリームの一番最初に、その構造を保つような変換を残すという処理をしてから後続処理を行う、という風に書いていました。でも、さらによく考えてみると、単結合のみの構造の列挙の段階で、似た処理を行っていることに気づきました。そこで、構造だけでなくその対称性も保存するような実装に変えて、その対称性を脱水素化ストリームでも流用するようにしました。これでトータル2秒ぐらい縮まりました。

flamegraph によるボトルネックの可視化

高速化のアイディアが出てこなかったので、ツールを使って分析することにしました。

FlameGraph と呼ばれる、プログラムの性能解析をしてくれるソフトウェアがあります。以下の紹介記事が簡潔かつ分かりやすいです。
qiita.com

こちらのツールは残念ながら Windows では使えないので、VirtualBox に Ubuntu を入れて、そこで実行しました。そのときに出てきたグラフがこんな感じです。

全体の横幅がプログラムの計算時間全体で、横幅の割合で各関数がそのうちどれぐらいの時間を使ったかが分かるようになっています。縦方向はコールスタックです。

この図からいろいろ読み取れますが、特に make_variants_symmetry という関数の実行時間がほぼほぼ全体の計算時間を決めていることが分かります。そのため、この関数の高速化に注力すればよいということが分かります。

比較の高速化

make_variants_symmetry を高速化すべく、先ほどの flamegraph を見ていると、SymmetricBitMatrix<N> 同士の比較にかなりの時間を使っていることに気づきました。これは内部的には u16 型が N 個（特に N=10）並んでいて、それらの比較を行っています。「そんなのが計算時間に効くのかよ」とビビりましたが、試しにこの部分を工夫してみました。

具体的には、比較のための情報圧縮を行いました。生データだと 16 bit×10 の比較になっていますが、実質は 10 × 10 の行列であり、しかも対角成分がすべて0の対称行列なので、その半分以下だけ比較できれば十分です。つまり実質 10 \times 9 \div 2 = 45 bit あれば十分なはずです。

そこで SymmetricBitMatrix<N> を64bit整数に変換してくれるメソッドを実装しました。一応、N \geqq 12 のときは64bitでは足りないので、128bitに変換するメソッドも準備しました。matrix.rs 内に impl_from_matrix_into_hash というマクロを用意して、u64, u128 それぞれに対して実行するという書き方をしています。

macro_rules! impl_from_matrix_into_hash {
    ($num_type:ty) => {
        impl<const N: usize> From<SymmetricBitMatrix<N>> for $num_type {
            fn from(mat: SymmetricBitMatrix<N>) -> Self {
                let mut hash = 0;
                let mut shift = 0;
                let mut mask = 0;
                for (i, &row) in mat.rows().iter().enumerate() {
                    hash |= ((row & mask) as $num_type) << shift;
                    shift += i;
                    mask = mask << 1 | 1;
                }
                hash
            }
        }
    };
}

impl_from_matrix_into_hash!(u64); // N <= 11 の場合のみハッシュが重複しないことが保証される
impl_from_matrix_into_hash!(u128); // N <= 16 で常にハッシュが重複しないことが保証される

自作の構造体同士の等値性の比較は PartialEq および Eq というトレイトで実装されます。これまでは derive マクロを使って自動実装していましたが、それをやめて、先ほどの情報圧縮を利用して比較する実装を自前で与えるようにしました。

impl<const N: usize> PartialEq for SymmetricBitMatrix<N> {
    fn eq(&self, other: &Self) -> bool {
        if N <= 11 {
            return u64::from(*self) == u64::from(*other);
        } else {
            return u128::from(*self) == u128::from(*other);
        }
    }
}

impl<const N: usize> Eq for SymmetricBitMatrix<N> {}

さらに、flamegraph の図によると HashSet への insert 部分で時間がかかっているとのことだったので、make_variants_symmetry で使っていた FxHashSet<SymmetricBitMatrix<N>> の代わりに、自前で set-like な構造体を実装しました。それが search.rs 内にある MatrixVariantSet です。

#[derive(Debug, Clone, Default)]
struct MatrixVariantSet {
    set_u64: FxHashSet<u64>,
    set_u128: FxHashSet<u128>,
}

impl MatrixVariantSet {
    fn insert<const N: usize>(&mut self, matrix: SymmetricBitMatrix<N>) -> bool {
        if N <= 11 {
            self.set_u64.insert(matrix.into())
        } else {
            self.set_u128.insert(matrix.into())
        }
    }

    fn contains<const N: usize>(&self, matrix: SymmetricBitMatrix<N>) -> bool {
        if N <= 11 {
            self.set_u64.contains(&matrix.into())
        } else {
            self.set_u128.contains(&matrix.into())
        }
    }
}

これのおかげで Set への挿入のタイミングでも圧縮された情報同士で比較してもらえるようになりました。これらを適用した後に再度 flamegraph を実行すると、以下のようになりました。Matrix 同士の比較の部分がごっそり小さくなったことが分かりました。計算時間もさらに2秒ぐらい短くなりました。

ノード番号入れ替えの高速化

make_variants_symmetry の中でもう一つ幅を利かせているのが、create_rearranged という SymmetricBitMatrix<N> のメソッドです。これは以下のように実装されています。

pub fn create_rearranged(&self, row_order: &[usize]) -> Self {
    let mut rows_new = [0; N];
    for (row_new, i_old) in rows_new.iter_mut().zip(row_order.iter()) {
        let row_old = unsafe { self.rows.get_unchecked(*i_old) };
        for (j_new, j_old) in row_order.iter().enumerate() {
            *row_new |= (row_old >> j_old & 1) << j_new;
        }
    }
    Self { rows: rows_new }
}

着目すべきはこの部分です。

unsafe { self.rows.get_unchecked(*i_old) }

通常配列の要素にインデックスでアクセスするときは、配列の外側のインデックスになっていないかのチェックを行い、もし境界の外だった場合は panic (プログラムを強制終了) するようになっています。

この境界チェック機能は、unsafe { } で挟まれたところでだけ使える get_unchecked というメソッドを使えばスキップすることができます*5。このメソッドが呼ばれる回数が多く、全体の処理時間に大きな影響を与えているので、この僅かなスキップが計算時間に影響するという訳です。実際、手元でこれを通常のインデックスアクセスに書き換えてみたところ、トータルの実行時間が1.7秒ほど伸びました。

ただ、この高速化は flamegraph を書いた時点ですでに行っていて、それ以外に高速化できそうな方法も見つからなかったので、断念しました。

おわりに

めちゃんこ高速化して楽しかったです（こなみ）
炭化水素の列挙はひと段落してしまったので、また何か面白いネタを探しておきます。ではまた。

追記：まさかの続編です。
smooth-pudding.hatenablog.com

*1:正確には炭素2個のみからなる四重結合の分子が存在するようです。 二原子炭素（C2）の化学合成に成功！ – 明らかになった4つの結合とナノカーボンの起源 | academist Journal ただこの一種類のみなので、今回はいなかったことにします。

*2:炭素数が2 のときだけ自動で満たされません(C2 のせい)。私のコードでは 炭素数が2 のときだけ和が3を超えないという条件に変えることで制御しています。

*3:dehydrogenate は脱水素化するという意味の動詞です。

*4:前回の記事より時間が短縮しているのは、前回の記事のコードで行っていなかった高速化の工夫を施したからです。同じ工夫を前回の記事のコードに適用すれば、おそらくもっと短い時間で終わるはずです。

*5:もちろん、スキップした分だけ、安全性を捨てることになります。get_unchecked で領域外にアクセスしたときの挙動は未定義動作 (コンパイラ依存 and/or やってみないと分からない) になります。絶対に領域内のインデックスにアクセスすることが保証されている状況でのみ使える技ですね。