另一種解鋼筋切割組合問題(無界背包問題)的方法(改自 thinker 所提觀念)

#!/usr/bin/env python
# -*- coding: utf8 -*-
""" cut(10, [7, 5, 3, 2]) = cut(3, [7, 5, 3, 2]) + cut(10, [5, 3, 2])
"""
def _plus1(a):
    a[0] += 1
    return a


def _cache(my_function):
    CACHE = {}
    def inner_function(*args):
        key = str(args[:])
        if not CACHE.get(key, None):
            CACHE[key] = my_function(*args)
        return CACHE[key]

    return inner_function


@_cache
def cut(total, sizes):
    propers = tuple([sz for sz in sizes if sz <= total])
    if len(propers) == 1:
        return [[0,]*(len(sizes)-len(propers))+[total/propers[0],],]
    elif not propers:
        if total < 0:
            return []
        else:
            return [[0,]*len(sizes),]

    result1 = [_plus1(a) for a in cut(total-sizes[0], sizes)]
    result2 = [[0, ]+a for a in cut(total, sizes[1:])]
    return result1 + result2


if __name__ == '__main__':
    from time import time

    bar = 10
    sizes = [7, 5, 3, 2]

    t0 = time()
    answer = cut(bar, sizes)
    print time() - t0

    print 'count: %s'%len(answer)
    #for a in answer:
    #    print a, sum([j*sizes[i] for i, j in enumerate(a)])

其計算邏輯是 cut(10, [7, 5, 3, 2]) = cut(3, [7, 5, 3, 2]) + cut(10, [5, 3, 2]) ，這觀念從 Thinker 那邊學來的，它用在計算 cut(10, [7, 5, 3, 2]) 的組合數有多少時，非常非常快，不過，我在擴充至求組合數有那些時，這方法的速度就比不上 cut(10, [7, 5, 3, 2]) = (cut(10, [5, 3, 2]) 的結果，其元素 0 插入 0) + (cut(3, [5, 3, 2]) 的結果，其元素 0 插入 1) 了。

不曉得是不是我在運算過程中，那裡少了優化?

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

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

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

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

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

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

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

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

閱讀更多…

讓 google 幫你作圖表

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

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

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

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

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

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

原來這就是背包問題呀!

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

那是一個解決鋼筋裁切的廢料量最佳化所衍生的難題。也就是要列出一根原料鋼筋要切成工地用尺寸的組合問題，如：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 的話，就停止處理。

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

== 後記 ==

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

陣列特殊排序需求

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

動態規劃問題

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

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


我們可以依序購買 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

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

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

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

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

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

須作到的工作有：

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

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

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

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

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

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

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 原理全忘了。我還是得一次只作一件事吧!