Union-Find書き直し(C#)

前々回のエントリで、Union-Findの効率が少し悪いとコメントをいただいた。確かに、実務上の問題は小さそうだが、このデータの持ち方はイマイチかもしれない。

そこで、次のように書き直してみた。

• 各ノードにParentとRank情報を持つ
• 配列ではなくツリー構造にする
• それに伴い、上限設定が不要に
• Intではなくobjectを受け付ける

objectなので、型が異なるもの同士もグループ化できる。

var uf = new UnionFind();

Assert.IsFalse(uf.IsSameGroup(1, "1"));
Assert.IsFalse(uf.IsSameGroup("2", 2));
Assert.IsFalse(uf.IsSameGroup("two", 2));

uf.Unite(1, "1");
uf.Unite(2, "2");
uf.Unite(2, "two");

Assert.IsTrue(uf.IsSameGroup(1, "1"));
Assert.IsTrue(uf.IsSameGroup("2", 2));
Assert.IsTrue(uf.IsSameGroup("two", 2));   //2と"two"が同グループになる
Assert.IsTrue(uf.IsSameGroup("2", "two"));

ソースコードはこちら。たぶん、ならし計算量がO(α(n))、ワーストケースがO(n)だろうか。

※α(n)：アッカーマン関数の逆関数と蟻本に書いてるがよくわからない。O(log(n))より速いらしい

    public class UnionFind
    {
        private class Node
        {
            public Node Parent { get; set; }
            public int Rank { get; set; }
        }

        private Dictionary<object, Node> _dict = new Dictionary<object, Node>();

        private Node Root(object data)
        {
            if (!_dict.ContainsKey(data))
            {
                var node = new Node();
                _dict.Add(data, node);
                return node;
            }
            else
            {
                var node = _dict[data];
                return Find(node);
            }
        }

        private Node Find(Node node)
        {
            if (node.Parent == null) return node;
            return node.Parent = Find(node.Parent);
        }

        public void Unite(IComparable x, IComparable y)
        {
            var xRoot = Root(x);
            var yRoot = Root(y);
            if (xRoot == yRoot) return;

            if (xRoot.Rank < yRoot.Rank)
            {
                xRoot.Parent = yRoot;
            }
            else
            {
                yRoot.Parent = xRoot;
                if (xRoot.Rank == yRoot.Rank) xRoot.Rank++;
            }
        }

        public bool IsSameGroup(IComparable x, IComparable y)
        {
            return Root(x) == Root(y);
        }
    }

Topcoderサイトのリンク集

ひどいデザインで有名なTopcoderのサイトだが、本当に目的の記事にたどり着けないレベルになってしまった。

個人的によく使うところをメモしておく。

SRM Match Editorials

http://apps.topcoder.com/wiki/display/tc/Algorithm+Problem+Set+Analysis

アクティブなマラソンマッチ一覧

https://community.topcoder.com/longcontest/?module=ViewActiveContests

過去のマラソンマッチ一覧

https://community.topcoder.com/longcontest/stats/?module=MatchList

マラソン結果

https://community.topcoder.com/longcontest/stats/?module=ViewOverview

Arenaのダウンロードリンク

http://www.topcoder.com/contest/arena/ContestAppletProd.jnlp

C#実装（Kruskal法とUnionFind）

クラスカル法（連結グラフの最小全域木を求める）はアルゴリズムがとても分かりやすくて好ましい。参考は例によって蟻本。

・辺の集合から、最小コストのものを見つける

・この辺が、これまで取り出した辺と連結していなかったら使う。連結済なら破棄

→これを辺がなくなるまで繰り返すだけ

連結の有無はUnionFindを使って判定するが、C#にはUnionFindがないので自前で用意した。連結リストで表現できればスマートだし、上限の指定も不要になるのだが、書くのが面倒なので配列を使うことに。

    /// <summary>
    /// union find, positive int only
    /// </summary>
    public class UnionFind
    {
        private int[] data;

        public UnionFind(int size)
        {
            data = new int[size];
            for (int i = 0; i < size; i++) data[i] = -1;
        }

        public bool Unite(int x, int y)
        {
            x = Root(x);
            y = Root(y);

            if (x != y)
            {
                if (data[y] < data[x])
                {
                    var tmp = y;
                    y = x;
                    x = tmp;
                }
                data[x] += data[y];
                data[y] = x;
            }
            return x != y;
        }

        public bool IsSameGroup(int x, int y)
        {
            return Root(x) == Root(y);
        }

        private int Root(int x)
        {
            return data[x] < 0 ? x : data[x] = Root(data[x]);
        }
    }

これを使って実装する。最小コストを順次取り出したいので、辺情報の格納にはSortedSetを利用した。リストでもいいのかもしれないが、ソートがQuckSortなので最悪ケースがO(N^2)になってしまうと考えた。

SortedSetが平衡二分木なので、計算量は辺の追加が平均・最悪ともにO(LogN)、コスト算出がほぼO(N)のはず。（UnionFindはほぼ定数と考えて良いらしい）

    /// <summary>
    /// Vertices index should be positive number
    /// </summary>
    public class Kruskal
    {
        SortedSet<ComparablePair<int, Tuple<int, int>>> pathSet = new SortedSet<ComparablePair<int, Tuple<int, int>>>();

        private int _maxIndex = 1;
        public Kruskal(int maxIndex)
        {
            _maxIndex = maxIndex;
        }

        public void AddPath(int s, int t, int cost)
        {
            var tuple = new Tuple<int, int>(s, t);
            pathSet.Add(new ComparablePair<int, Tuple<int, int>>(cost, tuple));
        }

        public int GetCost()
        {
            int ret = 0;
            var unionFind = new UnionFind(_maxIndex + 1);

            foreach (var path in pathSet)
            {
                int cost = path.Item1;
                int vertex1 = path.Item2.Item1;
                int vertex2 = path.Item2.Item2;

                if (!unionFind.IsSameGroup(vertex1, vertex2))
                {
                    unionFind.Unite(vertex1, vertex2);
                    ret += cost;
                }
            }

            return ret;
        }
    }

木星をバスケットボールとしたときの星の大きさ

宇宙少年の息子のために、ボールやビー玉で太陽系を作ってあげた。そのときにいろいろ計算した結果を、捨てるのも何なので日記に書いておく。

まず太陽系の恒星と惑星。木星の直系を24cm（大人用のバスケットボールくらい）に合わせると、地球と火星がビー玉、天王星と海王星がテニスボールくらい、土星が子供用のサッカーボールくらいだった。水星は小さいビー玉で、太陽の直径が部屋の床から天井までくらい。

距離感も再現しようとしたのだが、下記のとおりたとえボール大まで縮小しても、100km～1000kmレベルで離れてしまうので無理だった。宇宙ってかなりスカスカなんだな・・・。

※息子情報によると、土星はでかいくせに水に浮くらしい 

これが準惑星の結果。どれもゴマ粒である。そのくせに距離はものすごい。こうみると、冥王星を惑星カテゴリに入れるのは無理があるのが良くわかる。

小惑星、彗星と、衛星もやってみた。ゴマ粒～顕微鏡レベル。距離は面倒くさくなったので省略。

 ※小惑星と彗星の直系はかなり概算（球体でないので）

息子にせがまれて、銀河や宇宙の大きさも無理やりkmであらわしてみた。いずれも推計が諸説あるようなので、そのなかから適当に取り出した数値である。

ここまでくると、木星をバスケットボール大で表そうが、結局実感できないのであまり意味がない。

なお、数字マニアでもある息子は超大歓喜であった。

元Excelはこれ

yambe2002/astral_body_sizes · GitHub

二次平面のベクトル演算(C#実装)

二次平面幾何の復習として、内積や外積、線の交点などを計算するコードを書いてみた。どこかで再利用するかもしれないのでクラス化してある。

参考にしたのは蟻本とこのサイト。 

まず、X・Y座標の保持と、座標用の演算オペレーターを用意する。

    /// <summary>
    /// Pont/Vector class for 2D 
    /// </summary>
    public class Pt
    {
        public double X;
        public double Y;

        public Pt(double x, double y)
        {
            X = x;
            Y = y;
        }

        public static Pt operator +(Pt v1, Pt v2)
        {
            return new Pt(DoubleUtil.Add(v1.X, v2.X), DoubleUtil.Add(v1.Y, v2.Y));
        }
        public static Pt operator -(Pt v1, Pt v2)
        {
            return new Pt(DoubleUtil.Add(v1.X, -v2.X), DoubleUtil.Add(v1.Y, -v2.Y));
        }
        public static Pt operator *(Pt v, double d)
        {
            return new Pt(v.X * d, v.Y * d);
        }
    }

※なおDoubleUtil.Add()は加算のDouble誤差対応版で、実装はこんな感じ。

    public class DoubleUtil
    {
        private static double EPS = 1e-10;

        public static double Add(double a, double b)
        {
            if (Math.Abs(a + b) < EPS * (Math.Abs(a) + Math.Abs(b))) return 0;
            return a + b;
        }

        public static bool Eq(double a, double b)
        {
            return Math.Abs(a - b) < 1e-9;
        }
    }

次に、各点を原点Oからのベクトルとしたときの内積(dot product)と外積(cross product)を求めるメソッドを作成する。

内積：a・b ＝ |a| b| cosθ ＝ ax bx + ay by

外積：a × b ＝ |a| |b| sinθ ＝ ax by - ay bx

        /// <summary>
        /// Dot product of two vectors: O -> this and O -> other
        ///  a dot b = |a| |b| cos(theta) = ax bx + ax by
        ///  zero if two vectors run orthogonally
        /// </summary>
        public double Dot(Pt other)
        {
            return DoubleUtil.Add(this.X * other.X, this.Y * other.Y);
        }

        /// <summary>
        /// Cross(det) product of two vectors: O -> this and O -> other
        ///  a x b = |a| |b| sin(theta) = ax by - ay bx
        ///  zero if two vectors run parallelly
        /// </summary>
        public double Cross(Pt other)
        {
            return DoubleUtil.Add(this.X * other.Y, -this.Y * other.X);
        }

次は点qが線分p1-p2上にあるか判定するメソッド。

①ベクトルq→p1とq→p2が平行である ＝ 外積が0になる

②ベクトルq→p1とq→p2が反対を向いている ＝ cosθが負になる ＝ 内積が負になる

③qがp1またはp2と一致 ＝ cosθが0になる ＝ 内積が0になる

の条件で判定できる。

        /// <summary>
        /// point q exists on line p1-p2?
        /// </summary>
        public static bool OnSeg(Pt p1, Pt p2, Pt q)
        {
            return (p1 - q).Cross(p2 - q) == 0 && (p1 - q).Dot(p2 - q) <= 0;
        }

最後に、2直線の交点を求めるメソッドを作成する。

直前p1-p2上の点xは、変数tを使ってp1 + t(p2 - p1)と表現できる。これが直線q1-q2にのるための条件は、ベクトルq1→xがq1-q2上にあればよいので、この外積が0になればよい。

(q2 - q1) × (x - q1) ＝ 0

これからtをもとめ、次の条件式が導き出される。

p1 × ( (q2 - q1) × (q1 - p1) / ( q2 - q1 ) × ( p2 - p1 ) ) ( p2 - p1 )

        /// <summary>
        /// crosssing point of line p1-p2 and q1-q2
        /// </summary>
        public static Pt Intersect(Pt p1, Pt p2, Pt q1, Pt q2)
        {
            return p1 + (p2 - p1) * ((q2 - q1).Cross(q1 - p1) / (q2 - q1).Cross(p2 - p1));
        }

ただし、これは点が線分上にある条件① しか使っていない。そのため交点の有無を判定するメソッドを別に用意する必要がある。

線分q1-q2とp1-p2が平行な場合とそれ以外に分けて、OnSeg()とIntersect()を使って、点の位置関係ごとに判定する。

        public static bool HasIntersect(Pt p1, Pt p2, Pt q1, Pt q2)
        {
            if ((p1 - q1).Cross(p2 - q2) == 0)
            {
                return OnSeg(p1, q1, p2) || OnSeg(p1, q1, q2) || OnSeg(p2, q2, p1) || OnSeg(p2, q2, q1);
            }
            else
            {
                var r = Intersect(p1, q1, p2, q2);
                return OnSeg(p1, q1, r) && OnSeg(p2, q2, r);
            }
        }

これが完成版。2点間の距離を返すメソッドなども追加した。

    /// <summary>
    /// Pont/Vector class for 2D 
    /// </summary>
    public class Pt
    {
        public double X;
        public double Y;

        public Pt(double x, double y)
        {
            X = x;
            Y = y;
        }

        public static Pt operator +(Pt v1, Pt v2)
        {
            return new Pt(DoubleUtil.Add(v1.X, v2.X), DoubleUtil.Add(v1.Y, v2.Y));
        }
        public static Pt operator -(Pt v1, Pt v2)
        {
            return new Pt(DoubleUtil.Add(v1.X, -v2.X), DoubleUtil.Add(v1.Y, -v2.Y));
        }
        public static Pt operator *(Pt v, double d)
        {
            return new Pt(v.X * d, v.Y * d);
        }

        /// <summary>
        /// Dot product of two vectors: O -> this and O -> other
        ///  a dot b = |a| |b| cos(theta) = ax bx + ax by
        ///  zero if two vectors run orthogonally
        /// </summary>
        public double Dot(Pt other)
        {
            return DoubleUtil.Add(this.X * other.X, this.Y * other.Y);
        }

        /// <summary>
        /// Cross(det) product of two vectors: O -> this and O -> other
        ///  a x b = |a| |b| sin(theta) = ax by - ay bx
        ///  zero if two vectors run parallelly
        /// </summary>
        public double Cross(Pt other)
        {
            return DoubleUtil.Add(this.X * other.Y, -this.Y * other.X);
        }

        /// <summary>
        /// point q exists on line p1-p2?
        /// </summary>
        public static bool OnSeg(Pt p1, Pt p2, Pt q)
        {
            return (p1 - q).Cross(p2 - q) == 0 && (p1 - q).Dot(p2 - q) <= 0;
        }

        /// <summary>
        /// crosssing point of line p1-p2 and q1-q2
        /// </summary>
        public static Pt Intersect(Pt p1, Pt p2, Pt q1, Pt q2)
        {
            return p1 + (p2 - p1) * ((q2 - q1).Cross(q1 - p1) / (q2 - q1).Cross(p2 - p1));
        }
        public static bool HasIntersect(Pt p1, Pt p2, Pt q1, Pt q2)
        {
            if ((p1 - q1).Cross(p2 - q2) == 0)
            {
                return OnSeg(p1, q1, p2) || OnSeg(p1, q1, q2) || OnSeg(p2, q2, p1) || OnSeg(p2, q2, q1);
            }
            else
            {
                var r = Intersect(p1, q1, p2, q2);
                return OnSeg(p1, q1, r) && OnSeg(p2, q2, r);
            }
        }

        public static bool operator ==(Pt x, Pt y)
        {
            return (DoubleUtil.Eq(x.X, y.X) && DoubleUtil.Eq(x.Y, y.Y));
        }

        public static bool operator !=(Pt x, Pt y)
        {
            return (!DoubleUtil.Eq(x.X, y.X) || !DoubleUtil.Eq(x.Y, y.Y));
        }
        public double Norm()
        {
            return Math.Sqrt(X * X + Y * Y);
        }
        public double Dist(Pt other)
        {
            return (this - other).Norm();
        }
        public static double Dist(Pt v1, Pt v2)
        {
            return v1.Dist(v2);
        }
    }

ダイクストラ法の実装(C#)

前エントリでPriorityQueueが用意できたので、最短経路を解くコードをダイクストラで書いてみた。PriorityQueueを二分ヒープで実装したため、計算コストはO((E+V)logV)になっているはず。

ダイクストラ法 - Wikipedia

なお、C#のTupleがIComparableではないので、簡単なPairクラスも作成している。

    /// <summary>
    /// Get min cost between two points
    /// </summary>
    public class Dijkstra
    {
        private int maxIndex = -1;
        private const int INF = Int32.MaxValue;

        private int[,] _edgeArray;

        /// <summary>
        /// Dijkstra, get min cost between two points
        /// should not contain negatie cost path
        /// </summary>
        /// <param name="size">max index of vertices</param>
        public Dijkstra(int size)
        {
            maxIndex = size + 1;
            _edgeArray = new int[maxIndex, maxIndex];

            for (int i = 0; i < _edgeArray.GetLength(0); i++)
            {
                for (int j = 0; j < _edgeArray.GetLength(1); j++)
                {
                    _edgeArray[i, j] = INF;
                    if (i == j) _edgeArray[i, j] = 0;
                }
            }
        }

        // Add a path(no direction) with its cost
        public void AddPath(int s, int t, int cost)
        {
            _edgeArray[s, t] = Math.Min(_edgeArray[s, t], cost);
            _edgeArray[t, s] = _edgeArray[s, t];
        }

        //Get the min cost between s and t
        //return Int32.MaxValue if no path
        public int GetMinCost(int s, int t)
        {
            int[] cost = new int[maxIndex];
            for (int i = 0; i < cost.Length; i++) cost[i] = INF;
            cost[s] = 0;

            var priorityQueue = new PriorityQueue<ComparablePair<int, int>>(maxIndex);
            priorityQueue.Push( new ComparablePair<int, int>(0, s) );

            while (priorityQueue.Count() > 0)
            {
                var costDestinationPair = priorityQueue.Pop();
                if (cost[costDestinationPair.Item2] < costDestinationPair.Item1) continue;

                for (int i = 0; i < maxIndex; i++)
                {
                    int newCostToi = _edgeArray[costDestinationPair.Item2, i] == INF ? INF :
                        costDestinationPair.Item1 + _edgeArray[costDestinationPair.Item2, i];
                    if (newCostToi < cost[i])
                    {
                        cost[i] = newCostToi;
                        priorityQueue.Push(new ComparablePair<int, int>(newCostToi, i));
                    }
                }
            }

            return cost[t];
        }
    }

    public class ComparablePair<T, U> : IComparable where T : IComparable<T>
    {
        public readonly T Item1;
        public readonly U Item2;

        public int CompareTo(object obj)
        {
            ComparablePair<T, U> castedObj = (ComparablePair<T, U>)obj;
            return this.Item1.CompareTo(castedObj.Item1);
        }

        public ComparablePair(T t, U u)
        {
            Item1 = t;
            Item2 = u;
        }
    }

C#でのPriority Queue（優先キュー）実装

.NETのC#にはPriorityQueueがなくて困るので自作してみた。

二分ヒープを使った実装で、ほとんど蟻本のソースコードそのまま。ただし、次の拡張を行っている。

• ジェネリック型に対応
• 昇順、降順の指定可
• カスタムIComparerを指定可
• Peek()、Count()、Clear()、ToArray()を追加

平衡二分木への操作で、計算量はPushがO(logN)、PeekがO(1)。Clear()が実装上の都合でO(logN)になっている。

Linqは対応していないので、使いたい場合はToArray()で配列化する。

Githubはこちら。

https://github.com/yambe2002/topcoder_template/blob/master/MyLibrary.cs

public class PriorityQueue<T> where T: IComparable
    {
        private IComparer<T> _comparer = null;
        private int _type = 0;

        private T[] _heap;
        private int _sz = 0;

        private int _count = 0;

        /// <summary>
        /// Priority Queue with custom comparer
        /// </summary>
        public PriorityQueue(int maxSize, IComparer<T> comparer)
        {
            _heap = new T[maxSize];
            _comparer = comparer;
        }

        /// <summary>
        /// Priority queue
        /// </summary>
        /// <param name="maxSize">max size</param>
        /// <param name="type">0: asc, 1:desc</param>
        public PriorityQueue(int maxSize, int type = 0)
        {
            _heap = new T[maxSize];
            _type = type;
        }

        private int Compare(T x, T y)
        {
            if (_comparer != null) return _comparer.Compare(x, y);
            return _type == 0 ? x.CompareTo(y) : y.CompareTo(x);
        }

        public void Push(T x)
        {
            _count++;

            //node number
            var i = _sz++;

            while (i > 0)
            {
                //parent node number
                var p = (i - 1) / 2;

                if (Compare(_heap[p], x) <= 0) break;

                _heap[i] = _heap[p];
                i = p;
            }

            _heap[i] = x;
        }

        public T Pop()
        {
            _count--;

            T ret = _heap[0];
            T x = _heap[--_sz];

            int i = 0;
            while (i * 2 + 1 < _sz)
            {
                //children
                int a = i * 2 + 1;
                int b = i * 2 + 2;

                if (b < _sz && Compare(_heap[b], _heap[a]) < 0) a = b;

                if (Compare(_heap[a], x) >= 0) break;

                _heap[i] = _heap[a];
                i = a;
            }

            _heap[i] = x;

            return ret;
        }

        public int Count()
        {
            return _count;
        }

        public T Peek()
        {
            return _heap[0];
        }

        public bool Contains(T x)
        {
            for (int i = 0; i < _sz; i++) if (x.Equals(_heap[i])) return true;
            return false;
        }

        public void Clear()
        {
            while (this.Count() > 0) this.Pop();
        }

        public IEnumerator<T> GetEnumerator()
        {
            var ret = new List<T>();

            while (this.Count() > 0)
            {
                ret.Add(this.Pop());
            }

            foreach (var r in ret)
            {
                this.Push(r);
                yield return r;
            }
        }

        public T[] ToArray()
        {
            T[] array = new T[_sz];
            int i = 0;

            foreach (var r in this)
            {
                array[i++] = r;
            }

            return array;
        }
    }