以 Python 做測試驅動開發 (Test-driven development) 練習 -- 測試一個多邊形是否近似於四角形 (二)

回顧一下前篇文章中提到測試驅動開發(TDD)的三步驟:

1. 寫好程式的介面，只要能夠通過編譯即可
2. 增加測試
3. 執行測試，如果有失敗的測試的話就修改程式，然後回到步驟2

回到我們目前的IsQuadrilateral函數，目前的實作非常簡單，只是檢查傳入的多邊形是否有四個頂點，若是的話就傳回真。現在要來處理其他狀況:長的像四角形的五角形(或更多角)。解法如前篇所說，我們先取多邊形的第一個點做為p0，然後找出最遠的點p1，接著找離p1最遠的點p0。我們把這件工作交給FindFarthestPoint，此函數接受一個點p和一個集合S，傳回在集合S中離p最遠的點。這時我們依照TDD的精神，加上個測試先。testFindFarthestPoint利用我們之前定義好的square，測試離左下角最遠的點是否為右上角。

# 要放到 class TestSequenceFunctions 裡面
def testFindFarthestPoint(self):
    self.assertTrue( numpy.array_equal(FindFarthestPoint(numpy.array([0,0]), self.square)[0], numpy.array([1,1])) )
    self.assertTrue( numpy.array_equal(FindFarthestPoint(numpy.array([0,2]), self.square)[0], numpy.array([1,0])) )

def FindFarthestPoint(p, polygon):
    x = polygon[0,]-p[0]
    y = polygon[1,]-p[1]
    idx = numpy.argmax(x**2 + y**2)
    return polygon[:,idx], idx

好了，雖然繞了些路，終於萬事俱備，可以來寫判斷三角形的函數了。由於我們已經把多邊形切成兩半，各自是兩群連續的點的集合。我們可以想像從p1p2的連線，如果途中有一個凸出的點，p1到p2就會是一個三角形(或者更多角)。我使用CountMiddleVertex來計算一群連續的點之中有幾個凸出的點。這會是一個遞迴函數，先找到最突出的點p3，然後用p3切割線段，然後分別呼叫CountMiddleVertex繼續計算。

# 要放到 class TestSequenceFunctions 裡面
def testCountVertex(self):
    self.assertEqual(CountMiddleVertex(self.square2),   1)
    self.assertEqual(CountMiddleVertex(self.triangle),  1)
    self.assertEqual(CountMiddleVertex(self.square),    2)
    self.assertEqual(CountMiddleVertex(self.pentagon1), 3)
    self.assertEqual(CountMiddleVertex(self.pentagon2), 3)
    self.assertEqual(CountMiddleVertex(self.hexagon),   4)

def CountMiddleVertex(pointset, threshold=0.03):
    pointnum = pointset.shape[1]
    if pointnum <=2:
        return 0
    p1 = pointset[:, 0]
    p2 = pointset[:, -1]
    dist2 = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) 
    triarea = []
    for p in range(pointset.shape[1]):
        triarea.append(PolygonArea(numpy.column_stack((p1,p2,pointset[:,p]))))
    # 如果最突出的點超過門檻值
    if max(triarea)*1.0/dist2 > threshold:
        midpoint = triarea.index(max(triarea))
        vertexnum = 1 + CountMiddleVertex(pointset[:, 0:midpoint+1])
            + CountMiddleVertex(pointset[:, midpoint:], maxnum-vertexnum);
    return 0

CountMiddleVertex測試完成後，我們就利用它來判定三角形吧。更新後的IsQuadrilateral以及完整程式如下。為了讓CountMiddleVertex執行快一些，我加了一個參數maxnum，讓它可以早點結束。

import numpy
import unittest

def PolygonArea(polygon):
    x,y = numpy.column_stack((polygon, numpy.array([polygon[0,0],polygon[1,0]])))
    return abs((numpy.dot(x[0:-1],y[1:]) - numpy.dot(x[1:],y[0:-1])))/2.0 
    
def FindFarthestPoint(p, polygon):
    x = polygon[0,]-p[0]
    y = polygon[1,]-p[1]
    idx = numpy.argmax(x**2 + y**2)
    return polygon[:,idx], idx

def CountMiddleVertex(pointset, maxnum=10, threshold=0.03):
    pointnum = pointset.shape[1]
    if pointnum <=2:
        return 0
    p1 = pointset[:, 0]
    p2 = pointset[:, -1]
    dist2 = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) 
    triarea = []
    for p in range(pointset.shape[1]):
        triarea.append(PolygonArea(numpy.column_stack((p1,p2,pointset[:,p]))))
    if max(triarea)*1.0/dist2 > threshold:
        midpoint = triarea.index(max(triarea))
        vertexnum = 1 + CountMiddleVertex(pointset[:, 0:midpoint+1]);
        if vertexnum >= maxnum:
            return maxnum
        return vertexnum + CountMiddleVertex(pointset[:, midpoint:], maxnum-vertexnum);
    return 0

def IsQuadrilateral(polygon):
    x,y = polygon
    
    if len(x) < 4 or len(x)!=len(y):
        return False
   
    p1, p1idx = FindFarthestPoint(polygon[:,0], polygon)
    p2, p2idx = FindFarthestPoint(p1, polygon)
    
    if p1idx > p2idx:
        p1idx,p2idx = p2idx,p1idx
        p1,p2 = p2,p1
    
    if CountMiddleVertex(polygon[:, p1idx:p2idx+1], 2)==1 and \
        CountMiddleVertex(numpy.column_stack((polygon[:, p2idx:], polygon[:, :p1idx+1])), 2)==1:
        return True
    else:
        return False
    
class TestSequenceFunctions(unittest.TestCase):
    def setUp(self):
        self.triangle   = numpy.array( [ (0,10,0),(0,0,10)  ] )
        self.square     = numpy.array( [ (0,1,1,0),(0,0,1,1)  ] )
        self.square2    = numpy.array( [ (0,0.5,1,0),(0,-0.01,0,1)  ] )
        self.pentagon1  = numpy.array( [ (0,1,1,0,-1),(0,0,1,1,0.5)  ] )
        self.pentagon2  = numpy.array( [ (0,1,1,0,-0.01),(0,0,1,1,0.5)  ] )
        self.pentagon3  = numpy.array( [ (0,1,1,0,-0.21),(0,0,1,1,0.5)  ] )
        self.hexagon    = numpy.array( [ (1,1,0,-1,-1,0),(1,-1,-2,-1,1,2)  ] )
        
    def testFindFarthestPoint(self):
        self.assertTrue( numpy.array_equal(FindFarthestPoint(numpy.array([0,0]), self.square)[0], numpy.array([1,1])) )
        self.assertTrue( numpy.array_equal(FindFarthestPoint(numpy.array([0,2]), self.square)[0], numpy.array([1,0])) )
        
    def testCountVertex(self):
        self.assertEqual(CountMiddleVertex(self.square2),   1)
        self.assertEqual(CountMiddleVertex(self.triangle),  1)
        self.assertEqual(CountMiddleVertex(self.square),    2)
        self.assertEqual(CountMiddleVertex(self.pentagon1), 3)
        self.assertEqual(CountMiddleVertex(self.pentagon2), 3)
        self.assertEqual(CountMiddleVertex(self.hexagon),   4)
    
    def testIsQuadrilateral(self):
        self.assertEqual(IsQuadrilateral(self.triangle),   False)
        self.assertEqual(IsQuadrilateral(self.square),     True)
        self.assertEqual(IsQuadrilateral(self.pentagon1),  False)
        self.assertEqual(IsQuadrilateral(self.pentagon2),  True)
        self.assertEqual(IsQuadrilateral(self.pentagon3),  False)
        self.assertEqual(IsQuadrilateral(self.hexagon),    False)

if __name__ == '__main__':
    unittest.main()

以 Python 做測試驅動開發 (Test-driven development) 練習 -- 測試一個多邊形是否近似於四角形 (一)

這次用Python來練習驅動程式開發。
TDD簡單的說只有三步
1. 寫好程式的介面，只要能夠通過編譯即可
2. 增加測試
3. 執行測試，如果有失敗的測試的話就修改程式，然後回到步驟2

我要寫的程式是輸入一個多邊形的頂點，檢查是否為近似四邊形。
這個問題可以很簡單，我們只需要檢查多邊形的頂點數目，
如果是四個點的話就是四邊形。
當然我希望這個程式的能力不只如此，不然就沒有練習的效果啦。
其實這個問題是從實際的問題演變過來的，實際上是應用在擴增實境中，
我們從攝影機抓取影像，然後找出影像中的正方形，做為一個參考點。
一個平片上的正方形從不同角度去看，會變成一個任意的四邊形，
我們從影像分割出許多多邊形，然後判斷像不像四邊形。然而由於影像處理必定會有誤差，
所以一個四邊形可能有一些凸出的部分，而會變成五角形或六角形。我們的程式須要能處理這種情況。
這也是採用TDD的理由，我們可以一步一步改進我們的程式。

開始囉
首先定義函數IsQuadrilateral，輸入為一個多邊形，輸出真假值。由於現在還沒有測試，
不假思索就可以得到一個最簡單的實作方法，一律傳回假就好。

def IsQuadrilateral(polygon):
  return False

接著是增加測試，利用Python的unittest模組，我們將測試程式包在TestSequenceFunctions物件中。

import unittest
import numpy
class TestSequenceFunctions(unittest.TestCase):
  def setUp(self):
    self.triangle   = numpy.array( [ (0,1,0),(0,0,1)  ] )
    self.square     = numpy.array( [ (0,1,1,0),(0,0,1,1)  ] )
  def testIsQuadrilateral(self):
    self.assertEqual(IsQuadrilateral(self.triangle),    False)
    self.assertEqual(IsQuadrilateral(self.square),      True)
if __name__ == '__main__':
    unittest.main()

setUp函數是用來建立測試用的資料，這邊我們建立了兩個二維陣列來表示三角形和正方形的頂點。
numpy是常用的數學函式庫，程式中的 [ (0,1,0),(0,0,1)  ] 表示一個二維陣列，第一列是三個頂點的x座標 0,1,0，
第二列三個頂點的y座標 (0,0,1)。square物件存放的是正方形，有四個頂點。

完成以後就是執行測試。結果在第二個測試失敗了，所以我們要回頭修改程式。第二版的程式如下:檢查傳入的多邊形是否含有四個頂點，
有的話就為真，否則為假。

def IsQuadrilateral(polygon):
    x,y = polygon
    if len(x) != 4 or len(x)!=len(y):
        return False
    return True

現在回頭再跑測試，就可以過關啦。當然我們不能停在這裡，要繼續加測試才行。

五角形1:
self.pentagon1 = numpy.array( [ (0,1,1,0,-1),(0,0,1,1,0.5)  ] )
五角形2:
self.pentagon2 = numpy.array( [ (0,1,1,0,-0.01),(0,0,1,1,0.5)  ] )
五角形3:
self.pentagon3 = numpy.array( [ (0,1,1,0,-0.21),(0,0,1,1,0.5)  ] )
此處我們用了三個五角形，長相如下:

其中第二個五角形其實是四角形，只是左邊有一點小突出而已，我們的新測試如下:

def testIsQuadrilateral(self):
    self.assertEqual(IsQuadrilateral(self.triangle),   False)
    self.assertEqual(IsQuadrilateral(self.square),     True)
    self.assertEqual(IsQuadrilateral(self.pentagon1),  False)
    self.assertEqual(IsQuadrilateral(self.pentagon2),  True)
    self.assertEqual(IsQuadrilateral(self.pentagon3),  False)

修改演算法，四邊形是兩個三角形所組成，所以我們先嘗試找到四邊形的兩個頂點。這裡我用了一個經驗法則，先任意選取polygon上的一個點p0，然後選取離p0最遠的點p1，接著選取離p1最遠的點p2，如此一來p1和p2就是近似四邊形上的兩個點。p1p2的連線把polygon分成兩半，我們接下來只需要判斷這兩半是否為三角形即可。請見下面的圖。

為了要找到離p0最遠的點，我需要再寫一個函數FindFarthestPoint，同時需要新的測試來驗證這新的函數。
今天就先寫到這裡吧。來回顧一下到目前為止我們的程式看起來如何:

import unittest
import numpy

# 主函數
def IsQuadrilateral(polygon):
    x,y = polygon
    if len(x) != 4 or len(x)!=len(y):
        return False
    return True

# 測試程式
class TestSequenceFunctions(unittest.TestCase):
    # 此函數在 testIsQuadrilateral 被執行前會被呼叫，以建立測試所需的物件
    # 想一想為何需要一個獨立的 setUp 函數呢?
    def setUp(self):
        self.triangle   = numpy.array( [ (0,10,0),(0,0,10)  ] )
        self.square     = numpy.array( [ (0,1,1,0),(0,0,1,1)  ] )
        self.pentagon1  = numpy.array( [ (0,1,1,0,-1),(0,0,1,1,0.5)  ] )
        self.pentagon2  = numpy.array( [ (0,1,1,0,-0.01),(0,0,1,1,0.5)  ] ) # 這個測試會失敗
        self.pentagon3  = numpy.array( [ (0,1,1,0,-0.21),(0,0,1,1,0.5)  ] )

    def testIsQuadrilateral(self):
      self.assertEqual(IsQuadrilateral(self.triangle),   False)
      self.assertEqual(IsQuadrilateral(self.square),     True)
      self.assertEqual(IsQuadrilateral(self.pentagon1),  False)
      self.assertEqual(IsQuadrilateral(self.pentagon2),  True)
      self.assertEqual(IsQuadrilateral(self.pentagon3),  False)

if __name__ == '__main__':
    unittest.main()  # 會呼叫TestSequenceFunctions裡所有的測試函數

shuffling / random permutation

又到了野人獻曝的時間了 :)

最近和前輩在聊面試有那些問題可以出
前輩提到一個經典的問題：如何隨機排序一組撲克牌?
換個說法，就是要產生一種52張牌的排列。

前輩提到這個問題可以用排序演算法來解決。
我上維基查了一下，果然在Knuth的大作 The art of computer programming 中有給出兩個演算法 [1]。

假設我們要排序的陣列是 src，大小為N，以 0-based 作為索引方式。
第一個演算法是這樣的:產生另外一個陣列 a，給 a 中的每一個元素一個亂數，接下來對 a 排序，排序的同時將 src 內元素依相同順序移動。排序時可以選用任意的演算法，因此可以在 O(N lg N) 時間內完成。

Knuth 書中的第二個方法更有效率，只需 O(N) 時間即可完成。演算法是: 從 i=N-1 開始依次往前走訪陣列的每個位置，對於每個 i ，產生一個亂數 0<=j<i ，並交換 src[i], src[j]。以 python 來寫的話會像這樣:

for i in range(N-1, 0, -1):
    j = random.randint(0,i)
    src[i], src[j] = src[j], src[i]

有看過這些演算法的朋友會覺得到目前為止都是老調重彈，沒錯，請在往下看。我認為這兩個演算法其實是有關連的，第二個方法可視為第一個方法的特例: 也就是使用 selection sort 來做排序。selection sort 的想法為，每次從未排序的數列中挑出最大(或最小)的一個，將其移到排序完的數列。用簡單的例子來表現:

假設我們要按遞增順序排序這四個數字: 3 2 4 1

第一步是找出最大的數字 4，並放到陣列尾端
3 2 (4) [1] => 3 2 1 | 4

第二步會找到 3，並放到 4 的前面
(3) 2 [1] | 4 => 1 2 | 3 4

接下來會找到 2，順序不變
1 (2) | 3 4 => 1 | 2 3 4

排序完成!

觀察每一步的過程，可分成兩階段，第一階段是找出未排序數列中最大的數字，第二階段是交換到排序數列的開頭。其中第一階段原本要花費 O(n) 的時間以找出最大的值所在的位置，但是既然在第一個演算法中，每個位置的值都是亂數，所以我們可以直接從中隨機挑一個位置當作是最大的數字，這樣有相同的效果，而且時間複雜度只有 O(1)。

好，獻曝完畢 :p

以 selection sort 能夠說明第二個演算法是第一個方法的特例，也許嘗試其他排序法會帶來不同的發現呢。

[1] http://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle#The_modern_algorithm

Kurukshetra Athena 來自印度的數學程式競賽

上週從 Project Euler 傳來一個數學程式競賽的消息。
競賽的網址如下:
http://www.athena.kurukshetra.org.in/

這個競賽的題目都是數學相關問題，包含動態規劃、數論和許多領域，厲害的人也許能夠不依賴電腦就算出答案哩。雖然今年比賽已經結束了，有興趣的人還是可以去看看題目，並期待明年的比賽!

我今年沒有花很多時間，一周的時間內只完成了七題 :p
Score: 700 (完成七題)
Rank: 241

Question 6

There are 52735727364727372 students and 2889221271829121 teachers in a school . On children's day each teacher brings x toffees with him/her . All the toffees are collected and distributed equally to all the students . At the end it was found that there was exactly one toffee remaining .
What is the least possible value of x ?


我的想法
恩，這是很經典的一個問題。
首先我發現到分配之後只剩下一顆太妃糖，這表示老師的數目a和學生的數目b是互質的。接著把未知數 x 和 a,b 關聯起來，並且引入一個輔助變數 y 來表示每個學生所拿到的太妃糖數目，寫成數學形式是:
ax-by=1

接著我去查詢兩個互質的數有甚麼關連，發現了一個形式相似的公式: Bézout's identity，此公式說對於任意非零整數a和b，令其最大公因數為d，則存在x和y使得 ax+by=d。

要解出 x 和 y 可使用 Extended Euclidean algorithm，用 python 來寫的話只需要六行，而且還支援大數運算，真是相當輕鬆。

def extended_gcd(a,b):

  """Return x,y satisfy ax+by=gcd(a,b)"""
  if a % b == 0:
    return 0,1
  else:
    r = extended_gcd(b, a % b)
  return r[1], r[0] - (r[1] * (a / b))

找到 x, y 之後回到原來的問題。理論上符合 ax+by=d 的 (x,y) 有無限多組，此題目要求我們找的 x 為最小的正整數，剛好 extended_gcd 演算法得到的就是題目所要的。

加快遞迴函數的計算

演算法中有一個簡單但是實用的技巧叫 memoization，可用於加快遞迴函數的計算，基本的想法就是把計算過的值儲存起來，並在呼叫前檢查是否已經計算過。

比如說我們要計算Fibonacci數列的值，公式為

F(0) = 0

F(1) = 1

F(n) = F(n-1), F(n-2)

用C++來實作並配合memoization的話程式如下

int fibonacci(int n)
{
    if (n<2)
        return n;
    static map<int,int> mymap;
    map<int,int>::iterator it;
    it = mymap.find(n);
    if (it==mymap.end())
    {
        int f = fibonacci(n-1) + fibonacci(n-2);
        mymap[n] = f;
        return f;
    }
    else
        return it->second;
} 

主要的概念在於使用 map 來儲存計算過的 f(n)。

上述的方法需要自己在程式中檢查傳入的參數是否有計算過，
然而人是懶惰的，在某些程式語言中可以將一個函數直接"升級"，
獲得memoization的能力。
比如說下面是ypcat提供的python的範例

def memoize(f, cache={}):
  def g(*args):
    key = tuple(args)
    if key not in cache:
      cache[key] = f(*args)
    return cache[key]
  return g

def fib(n):
  return n if n<2 else fib(n-1)+fib(n-2)

fib = memoize(fib)
print fib(300) 

比C++更方便的地方在於使用了一個函數 memoize 來包裝遞迴函數，執行的時候還可以把 fib() 直接替換為包裝過的函數。

Python, C#, Java, Lisp 等支援 first-class function 的語言都能達成以上功能，
那C++有沒有辦法自動memoization呢? 我目前還不確定。