十月 1, 2011

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

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

之前的作法是將 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])

四月 27, 2011

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 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()

» 使用 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()

» 指派問題使用 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

八月 30, 2010

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 以 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.

八月 9, 2009
» The Fraction from a Decimal

定點數運算常用於 embedded systems 中,因為大部分低階的 MCU (例如: 8051, PIC, AVR 等)開發環境雖提供浮點運算,卻是軟體模擬的,除了慢,還明顯佔用原本就少得可憐的記憶體空間。 C/C++ 語言雖無定點數運算專用語法,程式員卻可通過手動調整,有效以整數運算完成相同效果。 定點數運作的原理,簡言之,就是將原來的實數(real number)或者小數(decimal),改寫成分數(fraction):如果 x 是個含小數的實數,我們可以找來兩個整數(p, q),將它們相除,來近似原來的 x (p/q ~= x)。 實務上人們可能還會要求上述的 q 要是 2 的冪次,因為電腦處理的都是 0, 1 的二進位運算, q 表示成 2 的冪次可以達到較高的精度;另一個原因我想是許多 f/w 程式員都患了 shift 偏執症 :p 有個友人前

七月 4, 2009
» The Analog Clock

……秒針急急忙忙的去撥動每一根短棒,使它們產生意義。然後分針慢吞吞的做同樣的事,使那些短棒產生另一種意義。三種針的位置和關係不斷變更,在錶面上切割出許多角來,夾住那不可捉摸的時間。……(摘自作文七巧:P86) 算一算日子,在現任公司混吃也有九個月了。很幸運的,一進來就參與一顆 ASIC 的開發,從一開始的寫 tools 測試 FPGA 功能,後來的寫 f/w 測試 ASIC ,到最後的參與產品開發。照規劃,一開始只打算拿來秀秀圖,偶爾也秀秀時間日期。後來為了把這顆小 MCU 的能耐完全壓榨出來,前些日子我還幫它加了類比鐘(Analog Clock)。自此,相框就不再僅僅只是相框了: 想起專科的畢業專題,我實作過一組函式庫,用來執行 3D 投影及相關的座標轉換。一晃眼已經十多年了,最近為了完成的這個類比鐘,竟然連描點畫線的程式都得自己手寫……

六月 16, 2008

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 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,這也就是我為什麼不用讀數學的關係,因為我還花滿多時間想數學的。

十月 10, 2007
» Not to be Fooled by Randomness

最近抱讀了一些機率的書,突然發現這門學問還滿有趣的。除了有趣,一些基本的機率分配及常見的隨機過程都非常實用。在複習的過程,當然也要作個紀錄: Random Number Generators 如果想跑些像 Monte Carlo 這類對隨機性很要求的 simulation ,有必要好好檢視一下手邊的隨機數產生器(Random Number Generator, RNG)是否合格。 線性同餘法(Linear congruential generators, LCGs)是多數開發環境標準的 RNG。它雖可以輕快、方便地產生隨機數;不幸地,它產生的隨機數,很容易就不夠隨機,不適合用在嚴肅的場合。 對於 RNG ,我們重視的有: 演算法要簡單、計算要快速 要能通過各種隨機性測試,特定模式不可出現在我們運用的場合。 用於密碼學領域時,還要夠安全,不易被破解

七月 20, 2007

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 數學證明:在相同個數的自然數子集合中,其子集合的積及和必屬惟一值

例如在 { 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 個數的自然數子集合推得。所以我們得知「在相同個數的自然數子集合中,其子集合的積及和必屬惟一值」。

biggo.com.tw

A Django site.