おねえさんを組み合わせ爆発から救う：完結編おねえさんは星になった

おねえさんを組み合わせ爆発から救うために、経路をZDDとして表したら、すっきりと経路情報が扱えました。
http://d.hatena.ne.jp/nowokay/20121018#1350528607

あとは、このZDDを効率よく構築できれば、おねえさんを救えそうです。このZDDの構築には、クヌース先生の開発したSimpathアルゴリズムを使うと非常に効率よく構築できます。

前回生成したZDDを見ると、同じノードにまとまっているものがいくつかあることがわかります。特に後半になるとどんどん同じパターンになるものがまとめられていきます。
つまり、この経路問題のZDDを構築するときには、いかに同じパターンになるものをまとめるかが鍵になるということです。

Simpathでは、辺の端だけに注目して、同じパターンになっていればそれ以降のノードを使いまわすという考え方で、ノードをまとめていきます。
つまり、この図では、1と8がつながっている、9と10がつながっている、という端のパターンが同じなので、これ以降のパターンは共通になると考えられます。

そこでSimpathでは、mate配列という接続情報保持配列を用意します。そして、ノードが端の場合は反対の端の番号、辺の途中の場合は0、辺が通らない場合はそのノード番号を保持します。

このmate配列のうち、すでに通っていて、それ以降も通る可能性がある部分をfrontierといいます。このfrontierをキーにして、一度計算したノードを保持しておくというわけです。

Simpathアルゴリズムについては、次のリンクからたどれるPDFにわかりやすく解説されています。
ZDDを用いたパスの列挙と索引生成
もともとは、The Art of Computer Programming 4 Fascicle1の演習問題225で解説されています。

The Art of Computer Programming Volume 4, Fascicle 1 Bitwise Tricks & Techniques; Binary Decision Diagarms 日本語版 (ASCII Addison Wesley Programming Se)
作者:Donald E. Knuth
発売日: 2011/05/20
メディア: 大型本

こうして実装すると、2x2マスの経路の場合、次のようなグラフが得られます。

ここで、前回のようなきれいなZDDではないことに注意してください。「6:8」「7:9」の節点は、既約のZDDではそれぞれ2個と1個ですが、ここでは3個と2個あります。
ZDDとしての処理が大切な場合は、後処理として、無駄なノードを削っていきます。この処理は、全節点数に比例した時間になります。今回は、経路数を数えるだけなので、そのまま扱います。

ところで、そう、このコードの実行時間です。
おねえさんのスーパーコンピュータで4時間半かかった8x8マスの経路数は6秒で、6年半かかった9x9マスの経路数は21秒で計算することができました！
そして、あの、おねえさんをロボットにしてしまった、25万年かけて計算した10x10マスの経路数は、なんと1分41秒で計算できました！290億年かかるといわれた11x11マスでも13分弱です。
これでおねえさんを救うことができそうです。

いや、しかし果たして、本当におねえさんを救えるのでしょうか？あの動画のおねえさんの決意を思い出してください。
『私は皆に「組み合わせ爆発のすごさ」を教えたいの！止めないで！』
高速なアルゴリズムを手にいれたら、おねえさんは、やはり限界まで経路を数えるんじゃないでしょうか。皆に「組み合わせ爆発のすごさ」を教えるために。
次のサイト*1によると、いま、この経路の問題の世界記録は、21x21マスのようです。やはりアルゴリズムにも限界はあります。
あの！お姉さんの努力を軽く凌駕・・湊離散構造処理系プロジェクトが道順の数え上げで世界記録を更新しました | キニナル！未来Lab．
結局、計算量は爆発してしまうのです。おねえさんのロボット化は止めれないように思います。

ところで、今回のプログラムで、生成されたグラフの節点数は、次のようになっています。

 1:8
 2:55
 3:303
 4:1471
 5:6646
 6:28753
 7:121102
 8:501318
 9:2051389
10:8327122
11:33608610

だいたい、4倍ずつ増えていっていることがわかります。そして単純に言えば、これは、使用メモリが4倍ずつになるということを表します。
今回のプログラムで、11x11マスを計算するとき、11GBのメモリを使いました。次のような実行オプションを指定しています。

> java -Xmx12288m -XX:NewRatio=5 simpath.Simpath

ぼくのマシンのメモリは16GBなので、これ以上の計算は無理です。時間の問題よりも、メモリの問題のほうがきつくなってくるのです。

さて、11x11マスのときで節点数はだいたい3x10^7になっています。また、4^10がだいたい100万なので、マスが10増えると、節点数は7桁増えることになります。ということは、60マス増えて71x71マスになったとき、節点数は42桁増えて、3x10^49くらいになると考えられます。そして地球の原子数は、だいたい3x10^50程度ではないかという話があります。原子ひとつでひとつのノードの情報を保持したとしても、71マスを計算するには地球と同じくらいの大きさのメモリが必要になってくるわけです。

しかし、おねえさんは、アルゴリズム以外の技術力は非常に高度なものがあります。自身をロボットにして25万年運用できるほどの！おねえさんなら、宇宙に飛び出して、星ほどの大きさのコンピュータを構築して、いつまでも計算し続けるのではないでしょうか。
もしかしたら、空に輝く一番明るい星は、経路数を計算するコンピュータになったおねえさん自身なのかもしれません。そのような星をずっと見ていると、たまにチカチカと激しくまたたくときがあります。おねえさんが、新しい計算結果をぼくらに教えてくれているに違いありません。

ありがとう、おねえさん、ぼくたちはおねえさんのことを忘れません。
おねえさん、永遠に・・・

まとめ
『フカシギの数え方』おねえさんといっしょにみんなも25万年数えよう
 おねえさんを組み合わせ爆発から救う：格子グラフを生成する
 おねえさんを組み合わせ爆発から救う：2分決定木を作る
 おねえさんを組み合わせ爆発から救う：無駄を省いてBDDを作る
 おねえさんを組み合わせ爆発から救う：ZDDがおねえさんを救う

package simpath;

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import kis.jgraphsample.CreateBdd.ThroughNode;
import kis.jgraphsample.CreateLattice.LatticeEdge;
import kis.jgraphsample.CreateLattice.LatticeNode;
import kis.jgraphsample.CreateTree.DiagramNode;
import kis.jgraphsample.CreateTree.EdgeNode;
import kis.jgraphsample.CreateTree.LeafNode;

public class Simpath {
    private static void connect(CreateTree.EdgeNode parent, DiagramNode child, boolean right){
        if(right){
            parent.right = child;
        }else{
            parent.left = child;
        }
    }
    
    public static void main(String[] args){
        long start = System.currentTimeMillis();
        //格子グラフの生成
        List<LatticeEdge> edges = new ArrayList<>();
        Map<Integer, LatticeNode> nodeMap = new LinkedHashMap<>();
        CreateLattice.generateGraph(12, nodeMap, edges);//数値は頂点数(マス数+1)
        
        int nodesize = nodeMap.size();
        //その辺以降で接続できるノードを得ておく
        Map<Integer, Set<Integer>> reachables = new HashMap<>();
        Set<Integer> reachable = new HashSet<>();
        for(int i = edges.size() - 1; i >= 0; --i){
            LatticeEdge edge = edges.get(i);
            reachable.add(edge.from.num);
            reachable.add(edge.to.num);
            reachables.put(i, new HashSet<>(reachable));
        }
        
        //mate配列準備
        int[] mate = new int[nodesize + 1];
        //初期化
        for(int i = 1; i < mate.length; ++i){
            mate[i] = i;
        }
        Map<String, DiagramNode> memo = Collections.synchronizedMap(new HashMap<String, DiagramNode>());
        
        List<List<DiagramNode>> layers = new ArrayList<>();
        for(int i = 0; i <= edges.size(); ++i){
            layers.add(Collections.synchronizedList(new ArrayList<DiagramNode>()));
        }
        int[] newmate = mate.clone();
        next(nodesize, edges, reachables, layers, 0, newmate, new HashSet<Integer>(), memo, true);
        
        System.out.println("count:" + count(root));
        int ns = 0;
        for(List<DiagramNode> l : layers){
            ns += l.size();
        }
        System.out.println("node:" + ns);
        System.out.println((System.currentTimeMillis() - start) / 1000);
        //CreateTree.paintDiagram(layers, 400, 300, 40, 50);
        //CreateTree.paintDiagram(layers, 450, 650, 40, 50);
        //CreateTree.paintDiagram(layers, 1200, 700, 35, 30);
    }

    static DiagramNode next(int nodesize, List<LatticeEdge> edges, Map<Integer, Set<Integer>> reachables, List<List<DiagramNode>> layers, int height, int[] mate, Set<Integer> reached, Map<String, DiagramNode> memo, boolean r) {
     
        DiagramNode left = new LeafNode(false);
        if(height == edges.size()) return left;

        List<DiagramNode> layer = layers.get(height);
        LatticeEdge edge = edges.get(height);
        Set<Integer> newreached = new HashSet<>(reached);
        newreached.add(edge.from.num);
        newreached.add(edge.to.num);
        
        //左側
        StringBuilder sb = new StringBuilder();
        if(!createFrontier(mate, reachables.get(height + 1), newreached, sb)){
            if(!memo.containsKey(sb.toString())){
                left = next(nodesize, edges, reachables, layers, height + 1, mate, newreached, memo, false);
                memo.put(sb.toString(), left);
            }else{
                left = memo.get(sb.toString());
            }
        }

        //右側
        int[] newmate = mate.clone();
        DiagramNode right = getRight(edge, nodesize, edges, reachables, layers, height, newmate, newreached, memo);
        
        //親ノードを作る
        boolean leftIsFalse = left instanceof LeafNode && !((LeafNode)left).value;
        if(right instanceof LeafNode && !((LeafNode)right).value){
            //右がfalse
            if(leftIsFalse){
                //左右ともfalseならfalseを返す
                return right;
            }
            ThroughNode tn = new CreateBdd.ThroughNode();
            connect(tn, left, r);
            layer.add(tn);
            return tn;
        }
        EdgeNode en = new EdgeNode(edge);
        if(right instanceof LeafNode && ((LeafNode)right).value){
            //右がtrue
            layers.get(height + 1).add(right);
        }
        en.right = right;
        if(!leftIsFalse){
            en.left = left;
        }
        layer.add(en);
        return en;
    }        
        
    private static DiagramNode getRight(LatticeEdge edge, int nodesize, List<LatticeEdge> edges, Map<Integer, Set<Integer>> reachables, List<List<DiagramNode>> layers, int height, int[] mate, Set<Integer> reached, Map<String, DiagramNode> memo) {
        
        int st = mate[edge.from.num];
        if(st == 0){
            //通過点にこれ以上つなげたらアウト
            LeafNode ln = new LeafNode(false);
            return ln;
        }
        int ed = mate[edge.to.num];
        if(ed == 0){
            //通過点にこれ以上つなげたらアウト
            LeafNode ln = new LeafNode(false);
            return ln;
        }
        
        if(st == edge.to.num && ed == edge.from.num){
            //ループができる
            LeafNode ln = new LeafNode(false);
            return ln;
        }
        
        //mate配列のつけかえ
        mate[st] = ed;
        mate[ed] = st;
        if(st != edge.from.num){
            mate[edge.from.num] = 0;
        }
        if(ed != edge.to.num){
            mate[edge.to.num] = 0;
        }

        if(mate[nodesize] == 1){
            //終点が始点とつながった
            LeafNode ln = new LeafNode(true);
            return ln;
        }
        //フロンティア接点から、キーを生成
        Set<Integer> reachable = reachables.get(height + 1);
        StringBuilder sb = new StringBuilder();
        if (createFrontier(mate, reachable, reached, sb)) {
            return new LeafNode(false);
        }
        //System.out.println(sb);
        DiagramNode next = memo.get(sb.toString());
        if(next == null){
            next = next(nodesize, edges, reachables, layers, height + 1, mate, reached, memo, true);
            memo.put(sb.toString(), next);
        }
        return next;
        
    }
    
    private static boolean createFrontier(int[] mate, Set<Integer> reachable, Set<Integer> reached, StringBuilder sb) {
        if(reachable == null) return false;
        for (int i = 2; i < mate.length; ++i) {
            if (reached.contains(i) && reachable.contains(i)) {
                if (mate[i] != 0 && //途中ではなく
                        mate[i] != 1 && //始点にもつながらず
                        !reachable.contains(mate[i])) //もうたどらないところにつながっていたらダメ
                {
                    return true;
                }
                sb.append(" ").append(i).append(":").append(mate[i]);
            }
        }
        return false;
    }
    
    /** カウント */
    static BigInteger count(EdgeNode node){
        return count(node, new HashMap<EdgeNode, BigInteger>());
    }
    static BigInteger count(EdgeNode node, Map<EdgeNode, BigInteger> cache){
        if(cache.containsKey(node)){
            return cache.get(node);
        }
        BigInteger c = BigInteger.ZERO;
        if(node.right != null){
            DiagramNode n = node.right;
            if(n instanceof LeafNode){
                return ((LeafNode)n).value ? BigInteger.ONE : BigInteger.ZERO;
            }
            c = c.add(count((EdgeNode)node.right, cache));
        }
        if(node.left != null){
            DiagramNode n = node.left;
            if(n instanceof LeafNode){
                return ((LeafNode)n).value ? BigInteger.ONE : BigInteger.ZERO;
            }
            c = c.add(count((EdgeNode)node.left, cache));
        }
        cache.put(node, c);
        return c;
    }
}

*1:検算もこのサイトの数字を使っています