PageRank



轉移公告

計劃把 http://blog.hoamon.info/ 文章全部轉移至 http://www.hoamon.info/blog/ 這裡,而本 Blogger 站台的文章近 500 篇,我預計在 2014-12-31 前移轉完畢,完成後 http://blog.hoamon.info/ 將只作代轉服務,一律把舊連結如 http://blog.hoamon.info/index.html 轉成 http://www.hoamon.info/blog/index.html ,敬請舊雨新知互相走告。

新文章只發佈在 http://www.hoamon.info/blog/

何岳峰 敬上
顯示具有 math 標籤的文章。 顯示所有文章
顯示具有 math 標籤的文章。 顯示所有文章

2011年9月16日 星期五

再改寫「背包問題」的求解程式碼

之前的作法是將 cut 函式所計算的 list 結果直接 append 到全域變數 tmps 中,這樣的 cut 函式是無法作 decorator 的。

新方法則是把 cut 函式的 input, output 重新規劃,讓答案就是 return 值,這樣 input 就能對應到單一 output ,透過這個特性,我們就能加上一個 @cache decorator ,去作快取。因為在求解的過程中,勢必會遇到重覆的 input ,有了快取,可以少算一次。

其中的 _no_cache_count 值指的是第一次遇到的 input 值,而 _cache_count 值則是利用 dictionary 找到答案的次數。

要怎麼建構出 cut 函式的樣貌? 我們一開始先抽象地想像這個 cut 函式要作到的事就是 answer = cut(bar, sizes)

answer 是我們要的答案,結構是 list of list( [[ , , , ..], [ , , , ..], ..] )。而 bar 是原長度, sizes 則是需求尺寸的 list 。

假設我們帶 bar = 10, sizes = [7, 5, 3, 2] ,那麼經過 cut 運算,就能得到一個 list of list 的 answer,那到底 answer 是多少? 我們先不管。但我們可以知道 10 拿給 7 去切,可以得到 0, 1 兩種組合。

所以 cut(10, [7, 5, 3, 2]) 一定會等於 cut(10, [5, 3, 2]) 的結果其解全部在元素 0 的位置插入 0 + cut(3, [5, 3, 2]) 的結果其解全部在元素 0 的位置插入 1 。

一樣地, cut(10, [5, 3, 2]) 也會等於 cut(10, [3, 2]) 的結果其解全部在元素 0 的位置插入 0 + cut(5, [3, 2]) 的結果其解全部在元素 0 的位置插入 1 + cut(0, [3, 2]) 的結果其解全部在元素 0 的位置插入 2 。

直到 cut(10, [2]) 時,我們知道它的結果就是 ([5], ) ,為什麼是一個 tuple of list ? 因為之前我們就定義 cut 一定要回傳 list of list ,而因為這次的 cut 回傳值本身並不會被修改,所以傳個 tuple 回去,可以少用一滴滴的記憶體(應該是一滴滴而已)。

當開始有 answer 被回傳後,我們就開始作合併的工作(就是把前一個需求尺寸的用量插入 answer 內的 list)。合併後再回傳。

程式碼如下,不過在實際跑的時候,有二件事我不能理解,為了比較 cut 與 cache_cut 的效率差別,我在同一個行程上分別跑了兩次 cut, 兩次 cache_cut ,而順序是 cache_cut, cut, cut, cache_cut ,cache_cut 比 cut 快,這很容易理解,但第二次的 cut 居然會比第一次的 cut 還慢,這我就不懂了。

另外,我每次跑 cut 之前,都是用 cs = CutSteel(bar, sizes) 創建新的物件,為什麼第二次跑的 cache_cut ,它還是可以找到第一次 cache_cut 所儲存的 CACHE 呢?

最後,我還得到一個結論,當解答組合數不多時,用 cut 會比 cache_cut 快。因為小題目,遇到重覆 input 少,但如果還是全部的 input 要儲存 CACHE ,那所花費的時間還不夠重覆 input 所節省的時間了。


1 #!/usr/bin/env python
2 # -*- coding: utf8 -*-
3 import types
4
5
6
7 class CutSteel:
8 u""" 目的:解鋼筋切割的組合問題(也就是背包問題),但不只是求組合數,
9 也要把所有的組合列出。
10 例: 10 公尺長的鋼筋,要切成 7, 5, 3, 2 公尺等,有多少種組合。
11 解:
12 7 5 3 2
13 [1, 0, 1, 0]
14 [1, 0, 0, 1]
15 [0, 2, 0, 0]
16 [0, 1, 1, 1]
17 [0, 1, 0, 2]
18 [0, 0, 3, 0]
19 [0, 0, 2, 2]
20 [0, 0, 1, 3]
21 [0, 0, 0, 5]
22 """
23 def __init__(self, bar, sizes):
24 if type(bar) != types.IntType or bar <= 0:
25 raise ValueError(u'只接受正整數')
26 for s in sizes:
27 if type(s) != types.IntType or s <= 0:
28 raise ValueError(u'只接受正整數')
29
30 self._no_cache_count = 0
31 self._cache_count = 0
32
33
34 def cache(my_function):
35 CACHE = {}
36 def inner_function(*args):
37 key = str(args[1:])
38
39 # try:
40 # #INFO 用 try 的會比 if 慢一點點。只慢一點點。
41 # CACHE[key]
42 # args[0]._cache_count += 1
43 # except KeyError:
44 # args[0]._no_cache_count += 1
45 # CACHE[key] = my_function(*args)
46
47 if not CACHE.get(key, None):
48 CACHE[key] = my_function(*args)
49 args[0]._no_cache_count += 1
50 else:
51 args[0]._cache_count += 1
52 return CACHE[key]
53
54 return inner_function
55
56
57 @cache
58 def bag(self, total, sizes):
59 u""" 只計算組合數 from thinker"""
60 propers = tuple([sz for sz in sizes if sz <= total])
61 if not propers:
62 if total >= self._minsize: return 0
63 else: return 1
64
65 num = self.bag(total - propers[0], propers) + self.bag(total, propers[1:])
66 return num
67
68
69 def cut(self, total, sizes):
70 u""" 本函式的 input 為「被切割長度」及「欲切割的種數」。
71
72 output 為該 input 的所有組合。
73 """
74 if len(sizes) == 1:
75 return (
76 [(total / sizes[0]), ],
77 )
78 elif total < sizes[-1]:
79 return (
80 [0,] * len(sizes),
81 )
82
83 return [
84 [j] + tr
85 for j in xrange(0, total / sizes[0] + 1)
86 for tr in self.cut(total - sizes[0] * j, sizes[1:])
87 ]
88
89
90 @cache
91 def cache_cut(self, total, sizes):
92 u""" 因為 cache_cut 函式本身是具有固定 input 就會產生固定 output ,
93 它們具有一對一或多對一的關係,所以我把 input,
94 output 放在一個 dictionary 中,若程式計算到相同的 input 時,
95 可免計算,直接從 dictionary 拿答案。
96
97 其實本函式就是複製 cut 函式後,
98 將函式內程式碼中的 self.cut 改成 self.cache_cut ,
99 並在函式名前加上 @cache 而已。
100 """
101 if len(sizes) == 1:
102 return (
103 [(total / sizes[0]), ],
104 )
105 # elif total < sizes[-1]:
106 # #INFO 多這個判斷式反而變慢了。因為已經用 cache 了,
107 # #所以那些 total < sizes[-1] 情況會變成比較少,
108 # #然而在一個 cache_cut 函式中多加一個 if ,則判斷時間會多一倍,
109 # #加速效果反而不如預期。
110 # return (
111 # [0,] * len(sizes),
112 # )
113
114 return [
115 [j] + tr
116 for j in xrange(0, total / sizes[0] + 1)
117 for tr in self.cache_cut(total - sizes[0] * j, sizes[1:])
118 ]
119
120
121
122 from time import time
123 import sys
124 if __name__ == '__main__':
125 #bar = sys.argv[1:]
126 #sizes = sys.argv[2:]
127 bar = 10
128 sizes = [7, 5, 3, 2]
129 sizes.sort(reverse=True)
130 sizes = tuple(sizes)
131
132 cs = CutSteel(bar, sizes)
133 cs._minsize = min(sizes)
134 print 'Total count: %s' % cs.bag(bar, tuple(sizes))
135
136 cs = CutSteel(bar, sizes)
137 time0 = time()
138 result = cs.cache_cut(bar, sizes)
139 print 'cache_cut spend time: %s' % (time() - time0)
140 print len(result)
141 print('\tno cache count: %s, cache count: %s'%(cs._no_cache_count, cs._cache_count))
142
143 # cs = CutSteel(bar, sizes)
144 # time0 = time()
145 # result = cs.cut(bar, sizes[:])
146 # print 'cut spend time: %s' % (time() - time0)
147 # print len(result)
148 # print('\tno cache count: %s, cache count: %s'%(cs._no_cache_count, cs._cache_count))
149 #
150 # cs = CutSteel(bar, sizes)
151 # time0 = time()
152 # result = cs.cut(bar, sizes[:])
153 # print 'cut spend time: %s' % (time() - time0)
154 # print len(result)
155 # print('\tno cache count: %s, cache count: %s'%(cs._no_cache_count, cs._cache_count))
156 #
157 # cs1 = CutSteel(bar, sizes)
158 # time0 = time()
159 # result = cs1.cache_cut(bar, sizes[:])
160 # print 'cache_cut spend time: %s' % (time() - time0)
161 # print len(result)
162 # print('\tno cache count: %s, cache count: %s'%(cs._no_cache_count, cs._cache_count))
163
164 for i in xrange(0, len(result)):
165 print(result[len(result)-i-1])

2011年4月13日 星期三

CMClass: 簡述 libsvm(Support Vector Machine library) 使用方法

libsvm乃台大林智仁老師開發的 Open source 工具,其目的為實作 Support Vector Machine 分類器,使用語言主要是 C++ ,目前也有 JAVA 版本,也提供其他語言的 wrapper ,像是 Perl, Python, Ruby, Matlab, Hashkell, Lisp 等。

詳細數學就不介紹了,怕大家睡著(但其實是因為還沒看懂),各位可以看一下下面那段這個影片,大略了解 SVM 分類器如何區別不同資料。



本文章主要介紹的是用 Python 語言去操作 libsvm 函式庫。

先解壓縮 libsvm.tgz 檔,可以看到 python 及 windows 資料夾,如果要在 Linux 中使用的話,請在主目錄中作

$ make lib

這樣會得到 libsvm.so.2 檔,這是 libsvm 的主函式庫,而在 windows 中使用的話,它則是先幫你編譯好這個檔了,可在 windows/ 找到這個 libsvm.dll 檔。

在 Linux 中,請把 python/*py 放到 /usr/local/lib/python2.6/site-packages 中,而 libsvm.so.2 放到 /usr/local/lib/python2.6/ 。

在 windows 中,請把 python/*py 放到 C:\Python26\Lib\site-packages 中,而 libsvm.dll 請放到 C:\Python26\Lib\windows 資料夾中(因為 svmutil.py 寫死了它的相對路徑,所以務必依它的相對位置置放)。

請在 Python shell 中,鍵入下列指令,測試是否安裝成功。

>>> from svmutil import *
>>>

沒錯誤訊息,那就是安裝對了。

使用 svm ,主要就是兩個動作: 訓練及預測。

訓練:

svmutil.svm_train 函式的引數有「類別標籤」、「觀察值」、「參數」。

你的原始資料若是如下:

1. 3, 4, 5, 6 => 第二類
2. 3, 4, 5, 5 => 第一類
3. ....

前面的 #. 表第幾個觀察值,後面逗號分隔的數據為各維度的值,行末則是放置該觀察值為第幾類的說明。請把它轉成

>>> Y = [2, 1, ...]
>>> x = [(3, 4, 5, 6), (3, 4, 5, 5), ...]

類別標籤請獨立放置到一個 list 中,而觀察值維度則依序放置到另一個 list 中。接下來,就能使用 svm_train:

>>> from svmutil import *
>>> model = svm_train(Y, x, '-c 4')

所得到的 model 就是一個經過訓練的分類器。

預測

接下來,我們要拿訓練好的分類器去預測新的觀察值:

>>> p_label, p_acc, p_val = svm_predict([0]*len(new_x), new_x, model)

而 p_label 就是依 new_x 順序所對應的類別標籤 list 。

下圖是我隨機生成的 300 點,圓點為原始的觀察值,而以線相連的連續點則是預測點。



詳細程式碼請參照如下:

 1 #! /usr/bin/python
2 # -*- coding: utf8 -*-
3
4 __author__="hoamon"
5 __date__ =u"$2011/4/12 下午 05:52:31$"
6
7 from math import pi, sin, cos
8 from random import random
9 from matplotlib import pyplot as plt
10 from svmutil import *
11
12 def circleData(centre, radius, down_limit_percent=0, lens=100, range=[0, 100]):
13 points = []
14 while len(points) < lens:
15 _angle = 2 * pi * random()
16 radius_percent = random()
17 if radius_percent < down_limit_percent: continue
18 _radius = radius * radius_percent
19 x = centre[0] + cos(_angle) * _radius
20 y = centre[1] + sin(_angle) * _radius
21 if range[0] <= x <= range[1] and range[0] <= y <= range[1]:
22 points.append((x, y))
23 return points
24
25
26 def test():
27 u""" 製作三群的隨機資料,每群皆 100 個點,點位置的 x, y 限制在 0 ~ 100 之間
28
29 最後利用 matplotlib 繪製出來的圖,"單點"表原始資料,而連續點畫線的部份,
30 該點位的類別則是利用 svm_predict 計算出來的。
31
32 Y = [1, 1, 1, ..., 2, 2, 2, ..., 3, 3, 3, ...]
33 x = [(x1, y1), (x2, y2), ...]
34 """
35 Y = [1] * 100 + [2] * 100 + [3] * 100
36 x1, x2, x3 = (circleData((35, 40), 12),
37 circleData((35, 40), 48, down_limit_percent=0.25),
38 circleData((80, 80), 20)
39 )
40 x = x1 + x2 + x3
41
42 m = svm_train(Y, x, '-c 4')
43
44 #INFO 在 100x100 的畫布上,打出 40000 個點,拿這 4 萬個點去給 m 作預測,算出這 4 萬個點的類別
45 points = [(i*0.5, j*0.5) for j in xrange(0, 200) for i in xrange(0, 200)]
46 p_label, p_acc, p_val = svm_predict([0]*40000, points, m)
47
48 line_1, line_2, line_3, pre_label = [], [], [], p_label[0]
49 for i in xrange(0, 200):
50 for j in xrange(0, 200):
51 index = i * 200 + j
52 now_label = p_label[index]
53 if now_label == 1 :
54 line_1.append(points[index])
55 elif now_label == 2 :
56 line_2.append(points[index])
57 elif now_label == 3 :
58 line_3.append(points[index])
59
60 fig = plt.figure()
61 ax = fig.add_subplot(111)
62 ax.plot([p[0] for p in x1], [p[1] for p in x1], 'ro')
63 ax.plot([p[0] for p in x2], [p[1] for p in x2], 'go')
64 ax.plot([p[0] for p in x3], [p[1] for p in x3], 'bo')
65 ax.plot([p[0] for p in line_1], [p[1] for p in line_1], 'r-', alpha=0.5)
66 ax.plot([p[0] for p in line_3], [p[1] for p in line_3], 'b-', alpha=0.5)
67 ax.set_title('Points of three classes')
68 ax.set_xlabel('x')
69 ax.set_ylabel('y')
70 ax.set_xlim(0, 100)
71 ax.set_ylim(0, 100)
72 plt.show()
73 return m, p_label, p_acc, p_val
74
75
76 if __name__ == "__main__":
77 test()

2011年1月23日 星期日

指派問題使用 python + lp_solve 解決

指派問題乃線性規劃的一種特例,它的特性是不須強調解為 0-1 變數或整數值,但最後算出來它卻一定就是 0 或 1 ,只是這種說法是學理上的,當我們使用程式來計算時,往往因為這些工具在計算過程中使用數值分析的方法,造成解的結果只會接近 0 或是 接近 1 ,而不是純正的 0 或 1。

這也就是我在第 73 行中,使用 v > 0 而不是 v == 1 的原因,如果是寫 v == 1 的話,有些 v 值是 0.999999 的,就不會顯現了。事實上,使用 v > 0.5 會更好。不過,我最後檢查時,發現 > 0 就可以秀出這 50 個 x 的值,也就算了。

lp_solve 函式庫安裝方法請見舊文


 1 # -*- coding: utf8 -*-
 2 """ 問題:
 3     
 4     指定 0, 1, 2, ..., 49 等 50 個不可重複的數字給 x0 ~ x49,例如 x0 = 12, x1 = 33, ...
 5
 6     y = sin(1*x0) + sin(2*x1) + sin(3*x2) + ... + sin(50*x49)
 7
 8     求解 y 之最大值?
 9
10     解法:
11
12     此問題可視為一種指派問題,也就是說有 0 ~ 49 等員工,要放到 x0 ~ x49 的職位去,
13     這樣決策變數就會變成 p00(值為 1 代表 x0=0), p01(值為 1 代表 x1=0),
14     p02 , ..., p49_49 等 2500 個決策變數,且其值必為 0 或 1 。
15
16     雖然目標函式看起來是非線性的,但其實是線性的, y 函式的係數應該長得如下:
17
18             x0          x1          x2          ...
19     0       0(C00)      0(C01)      0(C02)      ...
20     1       0.84(C10)   0.91(C11)   0.14(C12)   ...
21     2       0.91(C20)   -0.76(C21)  -0.28(C22)  ...
22     ...     ...         ...         ...         ...
23
24     所以如果決策變數是 p20 = p01 = p12 = 1,其餘為 0 ,則代表 x0 = 2,x1 = 0,x2 = 1,
25     這樣 y = 0.91 + 0 + 0.14 = 1.05 。
26
27     所以目標式可以寫成 y = C00 * p00 + C01 * p01 + ... + C49_49 * p49_49 。
28
29     最後再加上限制式
30     
31     p00 + p01 + ... + p0_49 = 1
32     p10 + p11 + ... + p1_49 = 1
33     ...
34     p49_0 + p49_1 + ... + p49_49 = 1
35
36     p00 + p10 + ... + p49_0 = 1
37     p01 + p11 + ... + p49_1 = 1
38     ...
39     p0_49 + p1_49 + ... + p49_49 = 1
40
41     等 100 條限制式後,就能求 y 的最佳解。
42
43 """
44 from math import sin
45 import lpsolve55 as L
46
47 LENGTH = 50
48 C = []
49
50 for i in xrange(LENGTH):
51     for j in xrange(LENGTH):
52         C.append(-1*sin((j+1)*i)) # lp_solve 預設解極小值問題,所以我把目標函數係數全乘以 -1
53
54 lp = L.lpsolve('make_lp', 0, LENGTH**2)
55 L.lpsolve('set_verbose', lp, L.IMPORTANT)
56 ret = L.lpsolve('set_obj_fn', lp, C)
57
58 for i in xrange(LENGTH):
59     p = [0,] * (LENGTH ** 2)
60     for j in xrange(i*LENGTH, i*LENGTH+LENGTH): p[j] = 1
61     ret = L.lpsolve('add_constraint', lp, p, L.EQ, 1)
62
63     p = [0,] * (LENGTH ** 2)
64     for j in xrange(0, LENGTH):
65         p[j*LENGTH+i] = 1
66     ret = L.lpsolve('add_constraint', lp, p, L.EQ, 1)
67
68 L.lpsolve('solve', lp)
69 print u'目標值: %s' % (L.lpsolve('get_objective', lp) * -1) #要乘以 -1 來還原目標值。
70 vars = L.lpsolve('get_variables', lp)[0]
71 print u'決策變數: %s' % vars
72 for (ij, v) in enumerate(vars):
73     if v > 0:
74         i = ij / LENGTH
75         j = ij % LENGTH
76         print 'x%s = %s, ' % (j, i),
77         if i % 5 + 1 == 5: print

目標值最佳解為 47.8620523191 。

各變數值如下:
x21 = 0, x32 = 1, x47 = 2, x33 = 3, x1 = 4,
x37 = 5, x16 = 6, x45 = 7, x11 = 8, x25 = 9,
x18 = 10, x30 = 11, x7 = 12, x17 = 13, x0 = 14,
x41 = 15, x36 = 16, x22 = 17, x49 = 18, x9 = 19,
x44 = 20, x26 = 21, x43 = 22, x13 = 23, x42 = 24,
x35 = 25, x8 = 26, x20 = 27, x39 = 28, x40 = 29,
x29 = 30, x10 = 31, x34 = 32, x4 = 33, x2 = 34,
x38 = 35, x24 = 36, x6 = 37, x46 = 38, x5 = 39,
x27 = 40, x28 = 41, x14 = 42, x23 = 43, x48 = 44,
x19 = 45, x31 = 46, x12 = 47, x15 = 48, x3 = 49,

=== 後記 ===

經老師指導後,使用

ret = L.lpsolve('set_binary', lp, [1,]*(LENGTH**2)) #大約加在第 59 行後

令決策變數為 0-1 二元變數後,計算時間馬上減少了 60% 。

道路施工的排程問題,類似背包問題,但須考慮不同的排列方式

如:長度 10 公尺的路面,若有 7 公尺/日、 5 公尺/日、 3 公尺/日的施工工班可供選擇,則有幾種的排程組合? 此問題很像背包問題,但在背包問題中,它不須考慮裝入物品的順序,而只考慮種類。若問題改為路面總長度 1000 公尺,而有 [260, 230, 190, 140, 80] 幾種工班時,我的解答是 69225 種。以下是我的解法:

 1 class LineSerial:
 2     u"""  目的:解路面排程問題,如:長度 10 公尺的路面,若有 7 公尺/日、 5 公尺/日、 3 公尺/日
 3         的施工工班可供選擇,則有幾種的排程組合。
 4
 5         解如下,共 17 種:
 6             [7, 7], [7, 5], [7, 3],
 7             [5, 7], [5, 5],
 8             [5, 3, 7], [5, 3, 5], [5, 3, 3],
 9             [3, 7],
10             [3, 5, 7], [3, 5, 5], [3, 5, 3],
11             [3, 3, 7], [3, 3, 5],
12             [3, 3, 3, 7], [3, 3, 3, 5], [3, 3, 3, 3],
13     """
14     def __init__(self, total, sizes):
15         """ serial_times 則是在計算 serial 函式被呼叫幾次。
16         """
17         sizes.sort(reverse=True)
18         self._sizes = sizes
19         self.serial_times = 0
20         self.result = []
21         self.serial(total, None, [])
22
23
24     def serial(self, total, length, tmp):
25         u""" 將 total 依序給 _sizes 中的所有元素去切,切完後就放入 tmp ,
26             當 total <= 0 時, 再放入 self.result 中。
27         """
28         #self.serial_times += 1
29         tmp.append(length)
30         if total <= 0:
31             self.result.append(tmp[1:])
32             return
33
34         for s in self._sizes: self.serial(total-s, s, tmp[:])
35
36
37
38 if __name__ == '__main__':
39     from time import time
40     total = 1000
41     lengths = [260, 230, 190, 140, 80]
42     time0 = time()
43     cs = LineSerial(total, lengths)
44     print time() - time0
45     print u'總組合數: %s' % len(cs.result)
46     print u'serial 遞迴次數: %s' % cs.serial_times
47     #for i in cs.result: print i

2011年1月22日 星期六

使用 genetic algorithm 來求解非線性問題。如 y = [ sin(1*x0) * sin(2*x1) ] + ... + [ sin(49*x48) * sin(50*x49) ],求 y 的最大值

問題描述: 指定 0,1,2,.........49 等50個數字給 x0~x49(不可重複),且

y = [ sin(1*x0) * sin(2*x1) ] + [ sin(3*x2) * sin(4*x3) ] + ... + [ sin(49*x48) * sin(50*x49) ]

請求解 y 之最大值?

本問題我使用 pyevolve 函式庫來幫我處理染色體的突變、重組、交配工作。我只需要提供目標函數(eval_func)的計算方式即可。本問題我的 y 最佳解是 20.4676 ,決策變數是 [33, 26, 36, 16, 45, 28, 37, 1, 19, 2, 25, 14, 0, 22, 6, 17, 35, 24, 11, 12, 27, 42, 49, 32, 13, 20, 23, 43, 41, 30, 4, 9, 21, 3, 10, 34, 38, 15, 18, 5, 47, 39, 44, 40, 8, 7, 31, 48, 46, 29] 。

 1 from pyevolve import G1DList
 2 from pyevolve import GSimpleGA
 3 from pyevolve import Selectors
 4 from pyevolve import Statistics
 5 from pyevolve import DBAdapters
 6 import pyevolve
 7 from math import sin
 8
 9 """ 指定 (0,1,2,.........49 等50個數字不可重複)給 x0~x49,例如  x0=12,  x1= 33, ....
10     y = [ sin(1*x0) * sin(2*x1) ] + [ sin(3*x2) * sin(4*x3) ] + ... + [ sin(49*x48) * sin(50*x49) ]
11     求解 y 之 最大值=?
12 """
13
14 def eval_func(chromosome):
15     score = 20.0 #為了不讓 score 的值小於 0,因為 pyevolve 不支援適存值小於 0 。
16     list = map(lambda a,b: (a, b), xrange(50), chromosome)
17     list.sort(key=lambda a: a[1])
18     for i in xrange(25):
19         score += sin((2*i+1)*list[i*2][0]) * sin((2*i+2)*list[i*2+1][0])
20
21     return score
22
23 # Enable the pyevolve logging system
24 pyevolve.logEnable()
25 # Genome instance, 1D List of 50 elements
26 genome = G1DList.G1DList(50)
27 # Sets the range max and min of the 1D List
28 genome.setParams(rangemin=0, rangemax=500)
29 # The evaluator function (evaluation function)
30 genome.evaluator.set(eval_func)
31 # Genetic Algorithm Instance
32 ga = GSimpleGA.GSimpleGA(genome)
33 # Set the Roulette Wheel selector method, the number of generations and
34 # the termination criteria
35 ga.selector.set(Selectors.GRouletteWheel)
36 ga.setGenerations(5000)
37 ga.terminationCriteria.set(GSimpleGA.ConvergenceCriteria)
38 # Sets the DB Adapter, the resetDB flag will make the Adapter recreate
39 # the database and erase all data every run, you should use this flag
40 # just in the first time, after the pyevolve.db was created, you can
41 # omit it.
42 sqlite_adapter = DBAdapters.DBSQLite(identify="ex1", resetDB=True)
43 ga.setDBAdapter(sqlite_adapter)
44 # Do the evolution, with stats dump
45 # frequency of 20 generations
46 ga.evolve(freq_stats=20)
47 # Best individual
48 print ga.bestIndividual()

2010年8月30日 星期一

以 least square method 搜尋非線性方程式的係數值

針對非線性迴歸求解其係數值,我寫了一份文件,考慮了其中繁雜的數學,我是用 rst 格式撰寫,配合 latex 轉成 PDF 。因為要轉 html 時,那些數學式實在是太麻煩了,所以我想就只提供 PDF, rst 給各位下載。

本文件用途說明:

當我們有一堆觀測數據,想要找出一條方程式可以將自變數、應變數的關係明確表示時,可以使用 least square method 求解,例如我們有 10 個觀察值(1, 4), (4, 5), (1, 2), (4, 5), (6, 9), (1, 8), (4, 5), (3, 4), (5, 5), (6, 6) 等,而我們先猜測它們的關係是二次的,應該要符合 Y = a * X^2 + b * X + c 方程式,有了觀測值及預估方程式形式後,接下來就可以利用 least square method 計算出 a, b, c 三個係數,期使它們能讓計算出來的 Y 值與觀察值的誤差最小。

在求取 a, b, c 三個係數上,我介紹了 4 種方法: the steepest descent method, Newton's method, the Guass-Newton method 和 levenberg-marquardt method 。

PDF表文下載

rst全文下載

** 以下純供搜尋引擎使用,生人勿近 **

應用問題說明
--------------------------------------------------------------------------------

當我們有一堆觀測數據,想要找出一條方程式可以將自變數、應變數的關係明確表示時,
可以使用 least square method 求解,例如我們有 10 個觀察值
(1, 4), (4, 5), (1, 2), (4, 5), (6, 9),
(1, 8), (4, 5), (3, 4), (5, 5), (6, 6) 等,
而我們先猜測它們的關係是二次的,應該要符合 :math:`$Y = a \times X^2 + b \times X + c$` 方程式,
有了觀測值及預估方程式形式後,接下來就可以利用 least square method 計
算出 :math:`$a$`, :math:`$b$`, :math:`$c$` 三個係數。

評估非線性方程式的係數值
--------------------------------------------------------------------------------

假定一含有 :math:`$n$` 個未知係數的
方程式 :math:`$E(x_1, x_2, ..., x_k, p_1, p_2, ..., p_n)$` ,
在 :math:`$m$` 組觀測值
(:math:`$(y_1, x_{11}, ..., x_{1k})$`, :math:`$(y_2, x_{21}, ..., x_{2k})$`, ..., :math:`$(y_m, x_{m1}, ...x_{mk})$`) 下,
欲求出此 :math:`$n$` 個係數的值,
期使 :math:`$E(x_1, x_2, ..., x_k, p_1, p_2, ..., p_n)$` 與觀測值的誤差最小。

首先定義理論值與觀測值的誤差如式 :ref:`\ref{errors define}` 的 :math:`$f_i$` 。

.. raw:: latex

\begin{equation}\label{errors define}
f_i(p_1, p_2, ..., p_n) = y_i - E(X_i, p_1, p_2, ..., p_n), i \in 1, ..., m
\end{equation}

\begin{equation}\label{F function}
F(p_1, p_2, ..., p_n) = \sum_{i=0}^{m} (f_i(p_1, p_2, ..., p_n))^2, m \geq n
\end{equation}

將每一觀測值誤差平方並加總起來,可得到一目標函式 :math:`$F$` 其定義如式 :ref:`\ref{F function}` 。
當 :math:`$ F $` 函式的值為全域最小值或是區域最小值,
其特定 :math:`$ [p_1^*, p_2^*, ..., p_n^*] $` 組合,
也就是 :math:`$P^*$` ,即為該 :math:`$ E $` 函式的係數全域或區域最佳解。

所有的非線性最佳化都是用迭代方法計算的:
從一個起點 :math:`$P_0$` 開始搜尋,產生一系列的向量 :math:`$P_1, P_2, ..., $` ,
希望可以收斂至一個 :math:`$ P^{*} $` ,使得函式成為全域或區域
最佳解 :cite:`\cite{k_madsen_imm_2004}` 。

.. raw:: latex

\begin{equation}\label{Taylor expansion}
F(P_{k+1}) = F(P_k) + h^{\top} {F}'(P_k) + \frac {1}{2} h^{\top} {F}''(P_k) h + O(||h||^3)
\end{equation}

\begin{equation}
h = P_{k+1} - P_k
\end{equation}

如式 :ref:`\ref{Taylor expansion}` ,首先以 Taylor expansion 改寫 :math:`$ F $` 函式,
:math:`$O(||h||^3)$` 代表後面無窮盡的項目,
:math:`$ k $` 代表的是第幾次 iteration 。在每次 iteration 中,
我們先找出 :math:`$ h $` ,使得 :math:`$F(P_{k+1}) < F(P_k)$` 。
目前常見計算 :math:`$h$` 的逼近法有 The Steepest Descent method,
Newton's method, Guass-Newton method, The Levenberg-Marquardt method。

The Steepest Descent method
********************************************************************************

Steepest Descent method 將式 :ref:`\ref{Taylor expansion}` 改寫
為式 :ref:`\ref{Taylor expansion2}` ,忽略二階以後項目。

.. raw:: latex

\begin{equation}\label{Taylor expansion2}
F(P_{k}+\alpha h) \simeq F(P_k) + \alpha h^{\top} {F}'(P_k), \alpha > 0
\end{equation}

而原來的 h 改寫成一個純量 :math:`$\alpha$` 再乘以向量,這樣在 :math:`$P_{k+1}$` 接
近 :math:`$P^{*}$` 時, :math:`$\alpha$` 必定近似 0 ,
所以我們又可以得到式 :ref:`\ref{steepest descent}` ,

.. raw:: latex

\begin{equation}\label{steepest descent}
\lim_{\alpha \to 0} \frac{F(P) - F(P + \alpha h)}{\alpha \left \| h \right \|}
= - \frac{1}{\alpha \left \| h \right \|} h^{\top} {F}'(P)
= - \frac{h^{\top}}{\alpha \left \| h \right \|} \frac{{F}'(P)}{\left \| {F}'(P) \right \|} \left \| {F}'(P) \right \|
= - \left \| {F}'(P) \right \| cos \theta
\end{equation}

:math:`$\alpha \left \| h \right \|$` 為純量,所以我們
希望 :math:`$F(P)-F(P+\alpha h)$` 的值愈大愈好,
這樣 :math:`$ cos \theta $` 就必須為 -1 ,
而 :math:`$cos \theta$` 是 :math:`$h$` 與 :math:`${F}'(P)$` 的夾角,
所以如式 :ref:`\ref{sd}` ,最佳的 :math:`$h_{sd}$` 必為 :math:`$-{F}'(P)$` 。

.. raw:: latex

\begin{equation}\label{sd}
h_{sd} = - {F}'(P)
\end{equation}

如式 :ref:`\ref{sd}` ,陡降法中的 :math:`$h_{sd}$` 是以斜率值負值為移動方向。
而 :math:`$\alpha$` 的值,
我們需用 line search 求得,但效率通常會是個問題,所以也可以使用 binary search 方式來求得,
其概念是先取得一個 :math:`$\alpha_{min}$` 值讓 :math:`$F(P)-F(P+\alpha h)$` 大於 0 ,
再取得一個 :math:`$\alpha_{max}$` 值讓 :math:`$F(P)-F(P+\alpha h)$` 小於 0 ,
接下來以 :math:`$\frac{1}{2}(\alpha_{min} + \alpha_{max})$` 為
新的 :math:`$\alpha_{middle}$` 值,
去計算 :math:`$F(P)-F(P+\alpha h)$` 是大於 0 還是小於 0 。若大於 0 ,則
新的 :math:`$\alpha$` 值
應為 :math:`$\frac{1}{2}(\alpha_{middle} + \alpha_{max})$` ,若小於 0 ,則新
的 :math:`$\alpha$` 應為 :math:`$\frac{1}{2}(\alpha_{middle} + \alpha_{min})$` ,
如此迭代計算後,當 :math:`$\alpha$` 滿足預設條件或達迭代次數即可決定。

Newton's method
********************************************************************************

牛頓法則考慮以 :math:`$ F $` 函式的二階 Hessian 矩陣來計算 h 。
它將式 :ref:`\ref{Taylor expansion2}` 再取其一次微分得到
式 :ref:`\ref{Taylor expansion for newton}` ,若 :math:`$(P_k + \alpha h) = P^{*}$` ,
則 :math:`${F}'(P_k + \alpha h )$` 必等於 0 。

.. raw:: latex

\begin{equation}\label{Taylor expansion for newton}
{F}'(P_{k}+\alpha h) \simeq {F}'(P_k) + \alpha h^{\top} {F}''(P_k), \alpha > 0
\end{equation}

所以我們可以得到式 :ref:`\ref{newton}` ,而 :math:`$h_{n}$` 就等於
式 :ref:`\ref{newton2}` 定義。

.. raw:: latex

\begin{equation}\label{newton}
0 = {F}'(P_k) + \alpha h^{\top} {F}''(P_k)
\end{equation}

.. raw:: latex

\begin{equation}\label{newton2}
h_{n} = - \frac{{F}'(P)}{\alpha {F}''(P)}
\end{equation}

在搜尋效率上,牛頓法為二元收斂,而陡降法為線性收斂,所以牛頓法在接近最佳解時比較快,
而陡降法則是離最佳解較遠時比較快,且因 Hessian matrix 在計算上
不一定為 positive definite ,所以牛頓法往往會混合陡降法來實作。

The Guass-Newton method
********************************************************************************

Least squares problems 一般能用陡降法或是牛頓法求解,但要追求效率的話,我們應該作部份調整。
像是盡量使用二元收斂或是不需實作出 :math:`$F$` 函式 :cite:`\cite{k_madsen_imm_2004}` 。

我們將式 :ref:`\ref{Taylor expansion}` 的 :math:`$F$` Taylor expansion 改使用
在 :math:`$f$` 上,如式 :ref:`\ref{f Taylor expansion}` 。

.. raw:: latex

\begin{equation}\label{f Taylor expansion}
f(x + h) = f(x) + J(x) h + O(||h||^2)
\end{equation}

\begin{equation}
(J(x))_{ij} = \frac{\partial f_i}{\partial x_j}(x) = {f}'(x)
\end{equation}

式 :ref:`\ref{f Taylor expansion}` 可再整理為式 :ref:`\ref{f Taylor expansion2}` 。

.. raw:: latex

\begin{equation}\label{f Taylor expansion2}
f(x + h) \simeq l(h) \equiv f(x) + J(x) h
\end{equation}

\begin{equation}\label{new F Taylor expansion}
F(x + h) \simeq L(h) \equiv \frac{1}{2}l(h)^{\top}l(h)
= \frac{1}{2}f^{\top}f + h^{\top}J^{\top}f + \frac{1}{2}h^{\top}J^{\top}Jh
\end{equation}

再將式 :ref:`\ref{new F Taylor expansion}` 對 h 作一次微分得到式 :ref:`\ref{gn function}` 。

.. raw:: latex

\begin{equation}\label{gn function}
{L}'(h) = J^{\top}f + J^{\top}Jh
\end{equation}

因為在最佳解時, :math:`${L}'(h)$` 等於 0 ,所以我們可以得到 :math:`$h_{gn}$` 如
式 :ref:`\ref{guass-newton}` 。

.. raw:: latex

\begin{equation}\label{guass-newton}
h_{gn} = - \frac{J^{\top}f}{J^{\top}J}
\end{equation}


The Levenberg–Marquardt method
********************************************************************************

The Newton's method may not be defined beacause of the singularity of
:math:`$ {F}''(P_k) $` at
a given point :math:`$ P_k $` , or the search direction :math:`$ h_n $`
may not be a descent direction
:cite:`\cite{bazaraa_nonlinear_2006}` .
The algorithm used the gradient methods their ability to converge from an
initial guess which may be outside
the region of convergence of other methods. The algorithm uses the Taylor
series method the ability to close
in on the converged values rapidly after the vicinity of the converged
values has been
reached. Thus, The method combines the best features of its predecessors
while avoiding their most serious
limitations. :cite:`\cite{marquardt_algorithm_1963}` .

Levenberg-Marquardt 方法乃加入一 damping parameter :math:`$ \mu $` 。

.. raw:: latex

\begin{equation}
h_{lm} = - \frac{J^{\top} f(P)}{J^{\top} J + \mu I}, J = {f}'(P), I = Identity\ Matrix
\end{equation}

此自定參數的效益在於,For all :math:`$ \mu > 0 $` the coefficient matrix is positive
definite,
and this ensures that :math:`$h_{lm}$` is a descent direction. For large
values of :math:`$ \mu $` we get

.. raw:: latex

\begin{equation}
h_{lm} \simeq - \frac{1}{\mu}{F}'(P)
\end{equation}

If :math:`$ \mu $` is very small, then :math:`$ h_{lm} \simeq h_{gn} $` ,
which is a good step in the
final stages of the iteration, when P is close to :math:`$ P^{*} $` .
If :math:`$ F(P^{*})=0 $` (or very small),
then we can get (almost) quadratic final convergence :cite:`\cite{k_madsen_imm_2004}` .

:math:`$\mu$` 值在每個 iteration 中,皆不同,而 :math:`$ \mu_{k} $` 的選擇,
主要有兩種方法,
一是看 :math:`$J(P_{k})^{\top} J(P_{k}) $` 中,
對角線元素中最大值再乘以 :math:`$ \gamma $` ,
一般來說, :math:`$\gamma$` 的值介於 :math:`$ 10^{-6} \sim 1 $` 之間。
或是也可以用 :math:`$ F(P_{k}) - F(P_{k-1}) $` 的值 :math:`$s_{k}$` 來判斷,
當 :math:`$ s_k \geq 0 $` 時
, :math:`$ \mu_k $` 增加 10 倍,
當 :math:`$ s_k \le 0 $` 時, :math:`$ \mu_k $` 減少 10 倍,
:math:`$ \mu_0 $` 的初始值則是設為 0.001 。

小結
--------------------------------------------------------------------------------

理論上,非線性最佳化是很難求得全域最佳解的,就算你運氣真的很好,
瞎貓怕上死耗子,讓你遇到全域最佳解,但你還是沒有辦法去驗證它的確是全域最佳解。
只有在某些特殊題型下,你才知道所求出的解是全域最佳解,最簡單的例子就一元二次方程式,
當它是凹向上時,有全域最小值,當它是凸向上時,有全域最大值。

所以,以上介紹的 4 種找尋 h 的方法,都只是一種有系統的找尋技巧,這些方法經過數學推導,
可算是比較有效率罷了,不代表只有這些方法可以用在 least square problem 上。
你也可以任意混合這 4 種方法,如在第 1 ~ 10 次迭代時,使用陡降法,在第 11 ~ 20 次時,
使用 Levenberg-Marquardt 方法,在第 20 次以後一律使用高斯-牛頓法,
這種方法是沒有對錯的,只有解題時間的快與慢而已,而這快與慢就又牽扯到起始值的挑選及非線性函式的特性。
這次這樣用比較快,不表示這對其他問題的解題速度就會比較快。

References
[1] M. S. Bazaraa, Hanif D. Sherali, and C. M. Shetty. Nonlinear programming. John Wiley
and Sons, 2006.
[2] K. Madsen, H. B. Nielsen, O. Tingleff, and Mathematical Modelling. IMM METHODS
FOR NON-LINEAR LEAST SQUARES PROBLEMS, 2004.
[3] Donald W. Marquardt. An algorithm for Least-Squares estimation of nonlinear parameters.
Journal of the Society for Industrial and Applied Mathematics, 11(2):431–441, June 1963.
ArticleType: primary_article / Full publication date: Jun., 1963 / Copyright 1963 Society
for Industrial and Applied Mathematics.

2009年10月4日 星期日

在工程專案中,我們如何評定資金的機會成本!

當我們評估一投資方案或是生產方案是否實行時,就現金流量而言,我們可用 IRR(Internal rate of return) 來作評估。然在比較機會成本上,試問你該以什麼樣的利率作為訂定報酬率的依據,選定 5% 而不是 3% 的理由是什麼呢?

用銀行報價的存款利率可以嗎? 那用那一家呢? 台新銀行、國泰銀行還是台中銀行? 試問萬泰銀行的存款利率高還是國泰銀行的利率高?

閱讀更多…

2008年6月15日 星期日

234+4378-2390+4012 = ?

遇到三位數以上的加減乘除時,我會直接把 python 直譯器打開,然後問它…

hoamon@ibmhoamon:~$ python
Python 2.5.1 (r251:54863, Oct 5 2007, 13:36:32) [GCC 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2)] on linux2 Type "help", "copyright", "credits" or "license" for more information.
>>> 234+4378-2390+4012
6234

我的腦袋已經不太對這類型的問題有耐心了。

小時候,在各學科的考試中,我只有數學是不須準備的,只要憑著老師上課教學的記憶及乖乖完成作業,我就可以上場考試了。而且在「數學」這個領域往往不會落後到三名之後。

還記得小一還是小二的時候,老師說要考九九乘法表,而在我還不知道什麼是九九乘法表時,已經有小朋友說,他幼稚園時已經背過了。記得那時候還覺得為什麼我的幼稚園跟別人不一樣,我只是在期待有熱狗、熱魚點心的下午趕快來到。

後來,老師解釋所謂的 2 x 2 就是把 2 加 2 次的意思,當我記住這個原則後,我就可以在腦中堆疊出 9 x 9 的答案,所以當我睡覺前還是沒事幹的時候(註1),我就開始回想九九乘法表,忘記答案的話,就用加法把它算出來。

記住原理比記住結果的樂趣還高,只是這種習慣讓我到了高中時,就遭受很大的挫折,記得是某次的高二物理考試,我花了一節課的時間,推導了一個公式,好讓我解一題 15 分的應用題,如此可想而知,我的高中成績應該是非常不理想的。

好了,要講我的重點了。每每看到了有關珠心算的新聞時,我就會想起九九乘法表的往事,到底我們的小朋友學會在腦中運行一個11位數以上加減乘除的目的在那裡? 打發時間、學好數學還是證明人比計算機強。

數學可以是很抽象地,就像埃米·諾特的對稱性定理,而這不是用珠算把所有 11 位數的數字作加減乘除練習後就會懂得。
  • 註1,這也就是我為什麼不用讀數學的關係,因為我還花滿多時間想數學的。

2007年12月29日 星期六

讓 google 幫你作圖表

過去作研究時,為了要看清楚數據的結果是不是與我們所想的一樣,我們會把數值轉成圖表的方式展示,這比起看到 1,2,4.0,8,... 的數字能更快速了解答案的正確與否。

因為以前都是用 matlab 作研究,所以生圖表時,是在單機上生成一個圖檔來看,若不用 matlab 也可以,把生成的數據倒給 gnuplot 一樣可行,或者是直接用 GNU Octave 也成。

我們現在都是用 python 了,而在 python 上也有十分優秀的圖表函式庫 matplot 可以用。

只不過,現在不只是要求數學程式化,我們也要作模式商業化,一個可行的解題方式,我們要讓使用者方便使用,最簡單的方法是讓它變成網站。這時候使用 matplot 就有點麻煩了。

如果可以使用 google chart api ,就會比較輕鬆,而且流量還可以丟過 google 處理。只不過,機密的數據還是不要透過 GET 方法讓 Proxy 儲存到,這時候,還是用 matplot 吧!

打個廣告,歡迎有程式設計能力且數學(離散、管理數學、統計)底子好的人才來我們 lab 唸碩博士,我們 lab 的目標主要是用資訊技術讓工程的品質、生產力提高。考試科目有營建管理概論、工程經濟/工程統計,基本上,我在土木系的大學部也只學過營建管理概論(三學分),工程經濟/工程統計是得自己另外唸的。

2007年12月20日 星期四

原來這就是背包問題呀!

日前(好久前了)提了一個讀研究所時所解決的問題,那時還滿自傲的,畢竟在讀土木的同學中,還沒看到有人會解。

那是一個解決鋼筋裁切的廢料量最佳化所衍生的難題。也就是要列出一根原料鋼筋要切成工地用尺寸的組合問題,如:18公尺的鋼筋若要切成 10, 7, 5, 4, 3 公尺長的鋼筋,則有那幾種切法。

Thinker 提了一個很好的方法來計算它的組合數有多少。不過,因為我要的是有那幾種組合,所以不曉得他的方法能否套用(當然可以)。

之前,我是用 Perl 來解這個問題的,不過,程式已經找不到了,不是我沒作版本控制,那時候用的是 CVS,然而在歷經多次系統安裝, CVS 的儲存庫已經不知在那了。

所以,我用 Python 重寫這個方法。有兩種解法,第一種是用我之前的觀念來解的。第二種是昨天想的,不過,第二種卻比較沒效率。
 1 #!/usr/bin/env python
 2 # -*- coding: utf8 -*-
 3 def cut(length, k, tmp):
 4     if k == sizeslen:
 5         if length >= sizes[-1]:
 6             return
 7         else:
 8             tmps.append(tmp[:])
 9             return
10
11     comp = int(length/sizes[k])
12     for i in xrange(comp+1):
13         j = comp - i
14         tmp[k] = j
15         cut(length-sizes[k]*j, k+1, tmp[:])
16
17 from time import time
18 import sys
19 if __name__ == '__main__':
20     bar = float(sys.argv[1])
21     sizes = [float(s) for s in sys.argv[2:]]
22     sizes.sort()
23     sizes.reverse()
24     sizeslen = len(sizes)
25     tmps = []
26     tmp = [0,]*sizeslen
27     time0 = time()
28     cut(bar, 0, tmp[:])
29     print time() - time0
30     print len(tmps)

第二種:
 1 #!/usr/bin/env python
 2 # -*- coding: utf8 -*-
 3 def cut(L, x, k, tmp, num):
 4     num[0] += 1
 5     diff = L-x[0]
 6     if diff > mat[0][0]:
 7         if tmp[x[1][0]] != 0: return
 8         else: tmp[x[1][0]] = x[1][1]
 9     elif diff < mat[0][0] and diff >= 0:
10         if tmp[x[1][0]] != 0: return
11         else: tmp[x[1][0]] = x[1][1]
12         tmps.append(tmp[:])
13         return 
14     else:
15         return
16     for (i, s) in enumerate(mat[k+1:]):
17         cut(L-x[0], s, i+k+1, tmp[:], num)
18         
19 def sort_by_value(k):
20     return (k, k[0])
21
22 from time import time
23 import sys
24 if __name__ == '__main__':
25     time0 = time()
26     bar = float(sys.argv[1])
27     sizes = [float(s) for s in sys.argv[2:]]
28     lensizes = len(sizes)
29     tmps = []
30     mat = []
31     for (i, s) in enumerate(sizes):
32         comp = int(bar/s)
33         for j in xrange(1, comp+1):
34             mat.append((s*j, (i, j)))
35     mat.sort(key=sort_by_value)
36     tmp = [0, ]*lensizes
37     num = [0]
38     for (i, m) in enumerate(mat):
39         cut(bar, m, i, tmp[:], num)
40
41     print time() - time0
42     print len(tmps)
43     print num[0]

新的觀念是把需求尺寸的倍數尺寸拿來當切割尺寸。如:10公尺要給7, 5, 3, 2 來切的話,我先把需求尺寸變成 7, 10, 5, 9, 6, 3, 10, 8, 6, 4, 2 等尺寸來切,如果在切的過程,剛好又遇到擁有相同因數的尺寸則跳過,像是10公尺已經被 8 切掉了,後來又遇到 2 的話,就停止處理。

在人工演練的過程中,覺得第二種方法所跑的迴圈數比較少,然而寫成程式後,效率卻比較差,且實際的迴圈數也比較多。

== 後記 ==

這裡有最新且更有效率的解法。

2007年9月28日 星期五

陣列特殊排序需求

Example: 排序 '140.120.109.1', '140.1.109.2', '14.120.109.1', '120.140.109.2'…等 ip 格式的文字串。這種排序的結果應是將 ip 切成 4 個部份後再來排序的。

>>> a = [4,5,6,1,2,3]
>>> a.sort()
>>> print a
[1,2,3,4,5,6]

上面的例子告訴我們,每個 list 都會有 sort 方法,而它的功能就是把陣列中元素作排序。如果陣列中元素有字串、另一個陣列、其他物件等, sort 方法一樣都可以排序,只是你對一個混合型態的陣列作排序,不容易有實質意義。

而如果陣列中的元素是單型態的,但它們的資料格式長得有點複雜,像是 ip 位置,那麼該如何排序呢!

>>> def sort_by_ip(X, Y):
... x_list = X.split('.')
... y_list = Y.split('.')
... for (i, x) in enumerate(x_list):
... if int(x) > int(y_list[i]):
... return 1
... else:
... return -1
... return 0
...
>>> ip = ['1.1.1.1', '140.120.109.1', '140.1.109.2', '14.120.109.1', '120.140.109.2']
>>> ip.sort(sort_by_ip)
>>> print ip
['1.1.1.1', '14.120.109.1', '120.140.109.2', '140.1.109.2', '140.120.109.1']

我 們利用一個自定的 sort_by_ip 方法來客製化我們的排序需求。陣列的 sort 方法,在執行時,會依序將陣列中的元素以兩個兩個為單位丟給 sort_by_ip 方法來處理,由 sort_by_ip 方法回傳出這兩個元素的大小,如果是前者為大,回傳值為 1 ,如果是後者為大,就回傳 -1 ,相等則回傳 0 。

所以你也就不須煩惱 sort 方法是不是用 quicksort/bubble 演算法,只須考量兩個元素的大小是如何判定的。

當問題只縮小到判定兩者大小時,就十分簡單了,在 sort_by_ip 方法中,我們首先拿到兩個元素的值,然後以 split 方法將其斷成 4 個子元素,再來依序比較這 4 個子元素的大小即可。

再來談一個例子:[{'cost': 10, 'id': 3}, {'cost': 1, 'id': 13}, {'cost': 100, 'id': 33}, {'cost': 84, 'id': 4}]該如何排序,解答如下:

def sort_by_cost(X, Y):
if X['cost'] > Y['cost']: return 1
elif X['cost'] <>
else: return 0

2007年9月18日 星期二

動態規劃問題

一個大問題本身是由許多小問題所組成,且求解小問題的概念,也就是求解大問題的概念,像這種大問題就可以使用動態規劃來求解,而事實上,我倒覺得動態規劃求解觀念就是遞迴觀念。

假設有一發電廠欲規劃三年內符合總需求機組數的每年購買機組數以使購買總成本最低,如下圖:


我們可以依序購買 3,1,3 ,其成本為6+1+10=17,也可以任序購買 3,2,2 ,其成本為3+4+5=12,如此我們可知 3,2,2 的購買方式比較好。而最佳解該如何求得?以最簡單想法來看,就是把每年的選擇窮舉出來即可,以上題計,只要找出 4^3 種方式,比較那些方式符合總需求機組數且成本最低即可。

但窮舉法總有一天跑不完,如果你的問題是有二十種選擇,且要規劃二十年的話,電腦就沒這能耐跑出來了。

所以此題的作法是把大問題切成小問題,我們不用管三年後該怎麼買,而是永遠只考慮下一年要怎麼買。

從第一年開始,因為至少購買一組,所以我們留下三種購買方式 {1, 2, 3} ,到第二年時,我們把第一年的 1 組,去加上第二年的可購買組數 [0, 1, 2, 3] ,可得到 {1, 2, 3, 4} ;把第一年的 2 組,去加上 [0, 1, 2, 3] 可得到 {2, 3, 4, 5};
依此觀念可得
{1, 2, 3, 4}
{2, 3, 4, 5}
{3, 4, 5, 6}
等方式,這也就是下一年(也就是第二年)的可能購買方式組合。接下來,我們把下一年(也就是第二年)的總數量為 X 的組合找出來,然後留下成本最低的那一種,結果我們可知下一年(也就是第二年)會留下 3, 4, 5, 6 等四種總數量的購買方式:

在第二年的總購買量為 3 時: 1, 2(成本為2+4); 2, 1(成本為3+1); 3, 0(成本為6+0) => 應選 2,1 這一組。
在第二年的總購買量為 4 時: 1, 3(2+5); 2, 2(3+4); 3, 1(6+1) => 剛好都一樣,隨便選 1, 3 吧!
在第二年的總購買量為 5 時: 2, 3(3+5); 3, 2(6+4) => 2, 3
在第二年的總購買量為 6 時: 3, 3(6+5) => 3, 3

所以我們得到在第二年總購買量可為 {3, 4, 5, 6} 而它們的購買方式為 (2, 1), (1, 3), (2, 3), (3, 3) 等四種。

再接下一年,我們把 {3, 4, 5, 6} 再分別與第三年的 [0, 1, 2, 3] 來作組合得到:

在第三年的總購買量為 7 時: (1, 3), 3(成本為7+10); (2, 3), 2(成本為8+5); (3, 3), 1(成本為11+2) => 所以選擇 2, 3, 2,其總成本為 13 。
在第三年的總購買量為 8 時: (2, 3), 3(8+10); (3, 3), 2(11+5) => 3, 3, 2(16)
在第三年的總購買量為 9 時: (3, 3), 3(11+10) => 3, 3, 3(21)

由此,我們選擇 2, 3, 2 為我們的三年購買方式,其總成本為最低,只有 13 。

以上是例題的簡單說明。

我們會把發電廠的完整問題,以 Python 程式解出。程式說明如下:

解答的資料結構:

Res = {
27: {
'cost': 40 , 'list': [5, 1, 3, 10, 2, 0, 0, 6, 0, 0],
'costlist': [10, 1, 4, 15, 2, 0, 0, 8, 0, 0]
},
....
}

27 代表當年度的總機組數量, cost 值代表在每年購買量為 [5, 1, 3, 10, 2, 0, 0, 6, 0, 0] 時的總成本,
costlist 陣列則是每年購買成本。

成本矩陣:

Cost = [
[0 , 2 , 5 , 7 , 9 , 10 , 12 , 16 , 18 , 19 , 20], # 第一年在各種購買量下的成本。
[0 , 1 , 5 , 7 , 9 , 12 , 15 , 18 , 21 , 25 , 27], # 第二年在各種購買量下的成本。
[0 , 2 , 3 , 4 , 7 , 9 , 10 , 13 , 15 , 19 , 22],
[0 , 1 , 2 , 6 , 7 , 8 , 10 , 11 , 12 , 14 , 15],
[0 , 1 , 2 , 4 , 5 , 7 , 11 , 12 , 14 , 16 , 18],
[0 , 2 , 3 , 6 , 9 , 11 , 12 , 13 , 16 , 17 , 21],
[0 , 2 , 6 , 7 , 8 , 9 , 11 , 14 , 16 , 20 , 22],
[0 , 2 , 4 , 5 , 6 , 7 , 8 , 12 , 14 , 17 , 19],
[0 , 1 , 3 , 5 , 8 , 9 , 10 , 12 , 13 , 17 , 18],
[0 , 3 , 7 , 10 , 12 , 15 , 18 , 21 , 24 , 25 , 27],
]

需求矩陣:

Limit = [3 , 4 , 9 , 11 , 14 , 18 , 18 , 21 , 24 , 27]

運算流程:

主要概念是把前一年度的購買組合與本年度的購買組作交叉組合配對,
會得到本年度的購買組合,每一購買組合的形式如下:

27: {
'cost': 40 , 'list': [5, 1, 3, 10, 2, 0, 0, 6, 0, 0],
'costlist': [10, 1, 4, 15, 2, 0, 0, 8, 0, 0]
},

27 代表當年度的總機組數量, cost 值代表在每年購買量為
[5, 1, 3, 10, 2, 0, 0, 6, 0, 0] 時的總成本,
costlist 陣列則是每年購買成本。

跑完當年度的購買組合後,就把不滿足需求量機組數的購買組合刪除。

等 10 個年度跑完後,再利用 sort_by_cost
函式作排序,把最便宜的購買組合放到最前面。

延伸閱讀:

sort_by_cost 函式的原理:請參考
http://wiki.python.org/moin/HowTo/Sorting#head-10e70070894a1bdb7580127b5cf764a44a2d6d29

for k_i in [k for k in Res.keys() if k < Limit[i]]: del Res[k_i] 的語法等價於:

toDel_k = []
for k in Res.keys():
if k < Limit[i]: toDel_k.append(k)

for k in toDel_k:
del Res[k]


完整程式下載: http://www.hoamon.info/_d/dynamic_programming.py

2007年7月20日 星期五

數學證明:S 為一有 k 個元素的自然數組合,其元素們的和、積依序為 a、b。若 T 亦為一有 k 個元素的自然數組合,其元素們的和、積亦依序為 a、b,則 S=T。

例如在 { 1, 2, 3, 4, 5, 6, 1, 1 } 自然數組合中,其積為 720,和為 23 ,是 { 1, 2, 3, 4, 5, 6, 1, 1 } 所獨有的,找不到其他同為 8 個數的組合,其積與和也等於 720 與 23 。

如果個數不同的話則有可能擁有相同的積與和,如: {3, 12} 與 {1, 1, 4, 9} 的積同為 36,和同為 15 。

證明方式如下:

首先證明在兩個個數的自然數組合中,其擁有獨特的積與和

在 X * Y = A ; X + Y = B 兩條方程式下,必定只有 1 或 2 個解或是沒有解。

本階段以函式圖說明。如下圖


兩線的交合點只有一個。而下圖


則有兩個交點,其解為 X=2, Y=3 或是 X=3, Y=2 ,但以組合的角度來看,這兩個解屬於相同的涵義。

基本上, X * Y = A 是一個凹向上的曲線,其於直線的交點只有 0, 1, 2 三種情形。又因為 X, Y 是我們任意選取的自然數,所以不會出現 0 個解的情形。所以我們可知在 2 個數的自然數組合中,必定擁有獨特的積與和。

接下來,證明在 3 個數的自然數組合中也具有相同的特性:

方程式改寫如下:
( X * Y ) * Z = A'
( X + Y ) + Z = B'
=>
A * Z = A'
B + Z = B'
=>
A * Z = A'
A + Z = B' + C
而 C = A - B

A 為自然數的積,所以它也會是自然數,於是我們得到 { A, Z } 的積 A' 與和 B' + C 是獨特的,除了 A 與 Z 之外,無法產生相同的積與和。

=> { X, Y, Z } 的積與和也會是獨特的。

最後 N+ 1 個數的自然數組合皆可由 N 個數的自然數組合推得。所以我們得知「相同個數的自然數組合,其和與積是惟一的」。

2007年3月8日 星期四

缺失報告及缺失編號的關係

課堂上,提到了一個工程案承辦人員時常面對的困擾:

當工程案被查核、督導後,就會有查核、督導委員的意見,這些意見有的簡短、有的複雜,所以承辮人員把這些意見分門別類到工程會制定的缺失編號時,往往傷透了腦筋,如何正確地快速地為這些意見作分類就是一個值得探討的問題。

初步想法是:
反轉一般搜尋的觀念,我用報告內容來找關鍵字。也就是說已經擁有一堆關鍵詞的集合,用這個關鍵詞集合來搜尋報告內容,然後依關連性高至低排列,讓使用者點選這篇文章應分類到那個集合。

系統並由此得到一次回饋,強化某報告與某集合的關連性。

須作到的工作有:
  1. 報告的中文斷詞
  2. 分析某分類的關鍵詞集合為何
  3. 關連的基本資料庫

自然數中的特定自然數集合?

大家似乎對數學問題的程式求解非常有興趣。我這邊再提供一個數學問題,這本來是我的碩士論文要作的題目(其中的一個部份)。

1條18公分長的竹筷,如果要分成10、7、2、3公分的話,且剩餘的長度不可大於等於2公分,這樣有那幾種分法?

程式要能寫到任意長度的竹筷可被任意個數的自然數切割。

我的程式目前只有 Perl 版的,過幾天改成 Python 版的,再 po 。

2007年3月3日 星期六

3X5 = 2X8 ??



上圖合理嗎?

答案請用滑鼠選取下面的空白處:

當然是不合理的。因為綠三角的斜率為 2/5 ,紅三角為 3/8 ,兩個合在一起怎麼會是1個直角三角形呢。

答案請用滑鼠選取上面的空白處。

微軟徵才試題

永源拿了兩個不同數字...
這兩個數字分別大於 1 ,也分別小於 50

永源只把這兩個數字的乘積告訴了亞譚 ...
永源再只把這兩個數字的和告訴了明歆 ...
永源問, 這兩個數字是什麼 ?

以下是亞譚和明歆的對話... (小強, 小毛, 小燕在旁坐陪)
亞譚: 我不知道這兩個數字是什麼!
明歆: 但我也還是不知道這兩個數字是什麼!
亞譚: 喔!那我知道那兩個數字是什麼了 !
明歆: 喔!那我也知道那兩個數字是什麼了 !
突然間 ....聰明的三位陪客同時也說: 我們也知道那兩個數字是什麼了 !
聰明的你, 告訴我....那兩個數字是什麼 ?

原文應是英文,翻成中文後,求解卻有多組,所以上文已有作過部份修改。

之前聽到這個題目時,是用質數的概念來解的,這方法不對。今天聽了老師的講法,才知道應是用集合的方式。

求解的邏輯如下:

找出兩個數字「積」的集合,把只出現1次的全部刪除,針對剩下數字對,把「和」的集合中,只出現1次的刪除。

這樣留下來的數字對只剩下638對。然後我們針對638對所產生的「積」,去找出對應數字對的「和」出現有幾次。如:
1344的積有兩個數字對(28,48)及(32,42),所以「和」為76及74,再去對應638對中的「和」是76的有幾次,74的有幾次。最後我們得到1344的「和」次數集合為(2,3)。

這一定不是我們要的答案,因為亞潭知道但明韶不知道的原因,在於只能有一個次數大於2的「和」次數集合,而上述的1344共有2個大於2的集合。

所以我們所要找尋的「和」次數集合應為(2,1,1,1,1)、(3,1,1,1)或是(6,1)之類的。

程式邏輯如下:

利用 dictionary(hash) 的特性,把「積」集合放在 times 變數,key為積值,value為對應的數字對陣列;「和」集合放在 plus 變數,key為和值,value為對應的數字對陣列。

刪完了只出現1次的數字對後,把剩下的「積」集合放到 s_times 變數;「和」集合放到 s_plus 變數中。

最後一個迴圈,則是把剩下的「積」集合中,所對應數字對的「和」次數集合(放到 a 陣列)找出來,並找出其次數只有1個大於2的 a 陣列,而 a陣列中惟一大於2的「和」其數字對即是解答。

Python 程式碼:

#!/usr/bin/python

(times, plus) = ({}, {})
for i in xrange(2, 50):
----for j in xrange(i + 1, 50):
--------if(times.has_key(i * j)):
------------times[i*j].append((i, j))
--------else:
------------times[i*j] = [(i, j)]
--------if(plus.has_key(i + j)):
------------plus[i+j].append((i, j))
--------else:
------------plus[i+j] = [(i, j)]
--------
(s_times, s_plus) = ({}, {})
for (t_k, t_v) in times.items():
----if(len(t_v) >= 2):
--------for (i, j) in t_v:
------------if(len(plus[i+j]) >= 2):
----------------if(s_times.has_key(i * j)):
--------------------s_times[i*j].append((i, j))
----------------else:
--------------------s_times[i*j] = [(i, j)]
----------------if(s_plus.has_key(i + j)):
--------------------s_plus[i+j].append((i, j))
----------------else:
--------------------s_plus[i+j] = [(i, j)]

for (t_k , t_v) in s_times.items():
----(a, b) = ([], 0)
----for (t_i, t_j) in t_v:
--------a.append(len(s_plus[t_i+t_j]))
----for a_i in a:
--------if(a_i >= 2): b += 1
----if(b == 1): print str(t_i) + ', ' + str(t_j)

後記:

老師講解這題的時候,是用 Matlab 程式求解的。我在課堂上,嘗試用 Python 求解,怎麼想都想不出來,直到回家路上,才發現看了 Matlab 程式後,一直把矩陣觀念放在腦中,而基本的 hash 原理全忘了。我還是得一次只作一件事吧!
Related Posts Plugin for WordPress, Blogger...