\( \newcommand{\ord}[1]{\mathcal{O}\left(#1\right)} \newcommand{\abs}[1]{\lvert #1 \rvert} \newcommand{\floor}[1]{\lfloor #1 \rfloor} \newcommand{\ceil}[1]{\lceil #1 \rceil} \newcommand{\opord}{\operatorname{\mathcal{O}}} \newcommand{\argmax}{\operatorname{arg\,max}} \newcommand{\str}[1]{\texttt{"#1"}} \)

2019年1月29日 星期二

[ Delaunay Triangulation ] Delaunay 三角剖分

我高中的時候在morris的網站第一次看到這個東西,那個時候就很想試著寫出來看看,但那時候英文很不好所以不會找資料就沒有成功。直到今年的NPSC,聽說台大出了一題找Voronoi Diagram的題目沒有人寫出來,於是我就想挑戰看看,寒假開始之後就認真研讀相關資料,終於在前幾天完成了隨機增量法以及分治法的模板。

先來說說甚麼是三角剖分 Triangulation 吧!假設平面上散布很多不重複的點,現在用一些線段將這些點連成許多的三角形,滿足所有三角形的內部不會有其他點,而且三角形的數量越多越好,這就是一個三角剖分。

Delaunay Triangulation是一種特殊的三角剖分,滿足所有的三角形其外接圓中不會出現其他的點,這樣的性質會使得所有三角形中最小的角盡量大,而且最重要的是它是 Voronoi Diagram 的對偶圖,Voronoi Diagram 直接求的話很容易有浮點誤差而且程式也不好寫,因此從Delaunay Triangulation直接轉換過去會更加方便。

還有一個很重要的點,透過Delaunay Triangulation 的邊集合可以找出平面歐幾里得距離最小生成樹。

這裡我先給出點的資料結構,方便接下來的講解:

一般來說找出 Delaunay Triangulation 的其中一個操作是要判斷一個點是不是在一個三角形的外接圓內,但直接去求外接圓圓心會有浮點誤差,因此有一個作法是將所有點投影到三維的拋物線曲面中,也就是下圖的碗狀曲面,會發現把一個圓投影上去之後,該圓會落在一個平面上,圓內的點會在平面下面,圓外的點會在平面上面,這樣就可以利用一些處理三維幾何的技術就可以避免浮點誤差的判斷點是否在外接圓內了。

範例程式碼中點在圓內輸出1,圓外輸出-1,圓上輸出0:

當然實作時並不會真正的寫一個三維點的資料結構,這樣程式碼會比較長。

Delaunay Triangulation 我所知道的有兩種,分別是隨機增量法以及分治法。隨機增量法再點隨機加入的時候複雜度是$\ord{n\log n}$,分治法則是穩定$\ord{n\log n}$,但我的實作分治法的效能是隨機增量法的三倍快。由於兩種做法都已經有詳細教學了這裡就不用文字敘述。

  1. 隨機增量法教學在《Computational Geometry - Algorithms and Applications》書中有提到,該書有中文版
  2. 分治法教學請點該網站連結
這裡附上程式碼,順帶一提morris的分治法好像在某些測資會因為三點共線爛掉

隨機增量法:

分治法:

2018年12月31日 星期一

[ Cyclic Longest Common Subsequence ] 環狀最長共同子序列

相信會看我發的文的人應該多少都知道LCS(最長共同子序列)是什麼吧?對於這樣的問題很早以前就有人提出$\ord{nm}$的方法了。
現在如果將問題中的字串改成頭尾相接環狀的字串,我們稱其最長共同子序列為CLCS,例如:
設A,B為頭尾相接環狀字串
A = "aabaa" , B = "aaabb"
CLCS(A,B) = "aaab"
像這樣的問題要怎麼解呢?剛好codeforces上面就有一題:http://codeforces.com/gym/100203/problem/B
為了解決這一題,查閱了不少奇怪的文章,最後給我找到了這篇論文:Solving Cyclic Longest Common Subsequence in Quadratic Time
對於CLCS提出了$\ord{nm}$的演算法,在這裡附上該演算法的實作:

作法大致上是將b字串變成原來的兩倍,然後做一般的LCS,但在過程中要注意更新答案的順序。接著修改DP來源表格,每次將來源表格的第一個row刪掉(reRoot函數),然後計算當前的答案。

當我AC這一題後,我看到了一個非常驚人的作法,也是用DP複雜度是$\ord{nm}$,但是我目前還無法完全理解為何這樣是可行的,而且這樣的做法我也不知道如何回溯找出一組解,只能知道CLCS的長度,希望路過的大大們能提供想法:

2018年12月23日 星期日

[ lower convex hull self-balancing binary search tree ] 下凸包平衡樹

下凸包平衡樹是一個動態維護下凸包的資料結構,他至少能在$\ord{\log n}$時間完成以下兩個操作:
  1. 查詢 point(x,y) 是否被包含於下凸包中
  2. 將 point(x,y) 加入至下凸包中,然後將不位於下凸包頂點上的點刪除
由於繼承自C++ STL map,因此所有能對map用的操作都能被使用,對於那些會改動map資料的操作請小心使用避免影響到凸包的正確性,例如operator[]就盡量不要使用。

對於某些動態規劃問題需要用凸包斜率優化時應該非常有用,有空會補上專門用來動態做斜率優化的版本

2018年11月30日 星期五

[ simplex algorithm ] 線性規劃 - 單純形法

線性規劃簡介

線性規劃是在線性約束條件下尋找某個目標線性函數的最大、最小值,例如:
$$
\begin{array}{lrl}
    \min& -3x_1+2x_2 & \\
    \mathrm{subject~to}
        &4x_1-6x_2&\le 5\\
        &5x_1-2x_2&\ge 7\\
        &x_1&\ge 0.
    \end{array}
$$

在最佳化問題中,線性規劃是非常重要的領域,很多的演算法問題都可以轉換成線性規劃問題來解決,例如:最大網路流、最大匹配等。因此研究如何系統化的解決線性規劃問題是非常重要的事情

線性規劃標準形

所有的線性規劃問題都可以轉換成下列標準型:
$$
\begin{array}{lrl}
    \max &\sum_{j=1}^n c_j \times x_j &\\
    \mathrm{subject~to}
        &\sum_{j=1}^n A_{i,j} \times x_j &\le b_j ,\; \forall i=1\sim m\\
        &x_j& \geq 0, \; \forall  j = 1\sim n
    \end{array}
$$

一般來說大多數解線性規劃的演算法都會要求將問題轉換成標準形,因此務必要理解其轉換方法。轉換成標準形的步驟如下:

  1. 將最小化改成最大化:
    只要改變目標函數的正負號,即可將最小化問題改成標準型設定的最大化問題,也就是說$$\min \sum_{j=1}^n c_j \times x_j$$等價於$$\max \sum_{j=1}^n -c_j \times x_j$$
  2. 去除等式:
    如果約束條件中存在等式,將其轉化為兩個分別為大於等於及小於等於的不等式取代,例如:
    $$x_1+x_2=5$$等價於$$x_1+x_2 \le 5\\x_1+x_2 \ge 5$$
  3. 去除大於等於:
    如果約束條件中存在大於等於約束,將約束兩邊取負即可,例如:
    $$x_1+x_2 \ge 5$$ 等價於 $$-x_1-x_2 \le -5$$
  4. 去除自由變數:
    如果未知數 $x_i$ 沒有非負約束,我們說 $x_i$ 是自由變數。加入新變量 $x_i'$,並用 $x_i-x_i'$替換原來的變量 $x_i$,並增加$x_i,x_i' \ge 0$的條件即可
標準化範例:
$$
\begin{array}{lrl}
    \min& -x_1+2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&= 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
首先將最小化改成最大化:
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&= 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
接著去除等式:
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&\le 5\\
        &x_1+x_2&\ge 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
去除大於等於:
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&\le 5\\
        &-x_1-x_2&\le -5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
最後去除自由變數,將$x_1$用$x_1-x_3$取代:
$$
\begin{array}{lrl}
    \max& x_1-2x_2-x_3 & \\
    \mathrm{subject~to}
        &x_1+x_2-x_3&\le 5\\
        &-x_1-x_2+x_3&\le -5\\
        &x_1-x_2-x_3&\le 2\\
        &x_1,x_2,x_3&\ge 0.
    \end{array}
$$

單純形法

輸入input

本人實作的單純形法其輸入是一個被稱為單純形表的 $m+1 \times n+1$ 矩陣$A$,標準形中每個項目的值都可以對應到矩陣$A$上:
$$
\begin{array}{lrl}
    \max &\sum_{j=1}^n A_{0,j} \times x_j &\\
    \mathrm{subject~to}
        &\sum_{j=1}^n A_{i,j} \times x_j &\le A_{i,0} ,\; \forall i=1\sim m\\
        &x_j& \geq 0, \; \forall  j = 1\sim n
    \end{array}
$$
可以發現:

  1. $A_{0,1} \sim A_{0,n}$對應標準形中的$c_1 \sim c_n$
  2. $A_{i,0}$對應標準形中的$b_i ,\; \forall i=1\sim m$
為了方便計算請設$A_{0,0}=0$

上面標準化範例寫成單純形表為以下形式:
$$
\begin{bmatrix}
0& 1 &-2 &-1 \\
5 &1& 1 &-1 \\
-5 &-1 &-1 &1 \\
2 &1 &-1& -1
\end{bmatrix}
$$

輸出output

輸出為一個$size=n+1$的vector $ans$,其中:
  1. $ans[0]$為線性規劃的解,也就是最大值
  2. $ans[1] \sim ans[n]$為變數$x_1 \sim x_n$的其中一組滿足條件的答案
  3. 如果無解或是解無限大則$ans$的$size=0$
  4. 注意如果答案是 $0$ 由於浮點數誤差的關係可能會產生 $``-0"$ 這種東西,要小心處理
上面標準化範例的解如下:
$$\{0.5,\; 3.5,\; 1.5,\; 0\}$$
意即在 $x_1=3.5,\;x_2=1.5,\;x_3=0$ 的情況下會有最大值 $x_1-2x_2-x_3=0.5$

程式碼

使用C++11

原理

有空時會補上,先讓我好好休息吧!終於可以正常睡覺了

2018年1月11日 星期四

[ Implement Queue using Stacks ] 用堆疊實作佇列

一般在C++ STL中queue都是用deque實作的,而deque用和vector類似的倍增法來實作(稍微多了一些操作導致常數有點大但是每次操作時間較為平均),雖然夠用但有些時候這不是一個好選擇。
大二資料結構期中考時,考了一題叫你用stack實作queue的方法,我那時候想了一下想了一個超帥的作法,需要用到兩個stack,我把它成為stack A和B。
  • push:
    push時直接把資料push進B中,$\ord{1}$
  • pop:
    這是個問題,我們把A的top當成是queue的front,所以直接把A pop就可以了,如果A是empty的話,就把B的資料依序pop然後push進A中,在進行該操作,均攤$\ord{1}$
  • front:
    如pop所述,把A的top當成是queue的front,如果A是empty再另做處理,$\ord{1}$
  • back:
    和front相反,把B的top當成是queue的back,如果B是empty再另做處理,$\ord{1}$
以下為實作過程,以c++ vector代替stack:

2017年10月18日 星期三

[ Median-of-Medians Algorithm ] 中位數演算法

這是一個可以在保證線性時間(c++ std::nth_element是隨機演算法)找出一個序列中第k大元素的演算法,網路上已經有不少教學,但是很多人都認為常數太大因此缺乏實作。

教學文在此:http://tmt514-blog.logdown.com/posts/484313-divide-and-conquer-method-iii-find-the-median

其實我高中時就想要試著去時做看看,但是因為那時的程式能力太差的關係,做出來的東西一直有bug,後來去忙其他事情後就被我忘掉了。最近因為學長面試有被問到一樣的問題跑來問我,才慢慢想起來有一份沒寫完的code,於是今天抱著不管怎樣都要寫出來的精神把他寫完了:


2017年10月13日 星期五

[ C++ std::sort python implement ] C++std::sort Python實作

最近有一個很靠北的課需要用python寫一些C++很容易做到的東西。
這次我們被要求寫一個quick sort,但是他要求的做法實在太糟糕了,於是我就參考了C++ std::sort來實做了一個quick sort(不包含heap sort的部分):