自己紹介
皆さんこんにちは。
山中に生息しているプログラマ＆SE、ほげふがーです。
計算量が普通なら爆発してしまうような数学の最適化問題を、量子コンピュータではなく、非力なPCやRaspbe〇〇piで、無理やり解かせるのが大好きです。
気になったら、何かコメントいただけると嬉しいです(^^
2020/05/04 msinit_shortest1.0.py…計算が実行されない問題を解決しました。
msinit_longest1.0.py …正しい最適解が計算されない問題を解決しました。
githubにアップしたもの
これです。
https://github.com/TiresiasGit/trunk/tree/master/mmiopt
git clone https://github.com/TiresiasGit/trunk.git
- randarray4.0.py
ex: python randarray4.0.py 10000x10000.npz 10000 uint8
In:ファイル名 Out:npzファイル(1~100の自然数の一様分布乱数の正方行列) - csv_to_npz.py
ex: python csv_to_npz.py InFile.csv OutFile.npz
─行列の対角要素は、全て0にして、csvカンマ区切りの表を作成してください。(excelなどで。)
一要素の入力可能な最大値は、9999999です。必要であれば、0~100の正規化をしてください。
In:csvファイル Out:npzファイル
下記のfilename.npzに、変換したOutFile.npzを入力してご使用ください。
-
msinit_shortest1.0.py
ex: python msinit_shortest1.0.py filename.npz txt 0.01 single no2-opt -
msinit_longest1.0.py
ex: python msinit_longest1.0.py 1.0.py filename.npz txt 0.01 single no2-opt
Input:filename int or float
Output:ansedge.npz or txt int
ansnode.npz or txt int
accumulation.npz or txt float
概要
最初は必ずMaximum Sum Initialize法で評価値(dataeva)を算出し,評価値が最大のノードをstart地点として決定します。
Maximum Sum Initialize法とMaximum Mean Initialize法との違い
Maximum Mean Initialize
(https://qiita.com/TiresiasGit/items/c6f2de7059e25ebf342a)
今までの自己最速記録に比べ、最大71倍高速化しました！
例：10000x10000のTSPにおいて
計算時間 50.9814秒→0.7162秒！
総距離 10417→10449 (わずかな悪化で済んでいる)
python msinit_shortest1.0.py 10000x10000.npz txt 0.01 single no2-opt
start load file!
[[ 0 68 82 ... 12 67 18]
[70 0 98 ... 14 36 85]
[53 99 0 ... 3 10 77]
...
[ 4 15 64 ... 0 77 80]
[80 33 64 ... 11 0 36]
[92 5 16 ... 1 80 0]]
end load file!
flags = [0 0 0 ... 0 0 0]
start max mean initalize.
er=0.01
n=3137
end max mean initialize.
start search min edges.
end search min edges.
start make ansedge from ansnode.
end make ansedge from ansnode.
search end ansnode-> [7133 35 66 ... 9976 9839 9728]
Accumulation Distance = 10449.0
ansedge = [ 1. 1. 1. ... 45. 79. 26.]
calculation elapsed time=0.7161514759063721
total elapsed time=2.8556344509124756
>python mmiopt_shortest3.0.py 10000x10000.npz txt
start load file!
end load file!
start max mean initalize.
end max mean initialize.
start 2-opt
end 2-opt
search end ansnode-> [8132 7417 211 ... 7404 2487 5113]
Accumulation Distance = 10417.0
ansedge = [ 1. 1. 1. ... 11. 81. 29.]
calculation elapsed time=50.98143482208252
total elapsed time=53.1019127368927
評価値の計算
Maximum Mean Initialize法 Maximum Sum Initialize法
・平均値 → 合計値（開発のなりゆきでこうなりました...）
・行の内の全要素がInput → 統計学的なサンプルサイズのみInput
とよりヒューリスティックス（直感型）に仕様変更しました。
最小値の探索において探索済み除外の改善
Maximum Sumに変更したのとは別で、探索済み除外の改善をしました。
従来 mmiopt_ver3.0
→リスト内包表記とflagsの組み合わせ
while num > 0:
：
indexes = [i for i, x in enumerate(data[index]) if 0 == flags[i]]
最新 msinit_ver1.0
→flagsの探索済み(index)に100を代入、np.add()でdata[index]に行単位で加算する(各ノードの行に対して一度だけ適用されるため数値的発散の心配はない)
while num > 0:
：
data[index]=np.add(data[index],flags)
:
flags[index] = 100 #訪問済みにする
10000x10000のTSPにおいて、
トータル計算時間50.7954秒→25.7188秒
となりました。TSPでは最小値のサーチ時間が計算時間の9割以上を占めるので,その中の探索済みにflagを立てる処理を、上記処理に変更したことにより、その部分が2倍速ぐらいとなり、全体としてほぼ2倍速となりました。
ハイパーパラメータの考察
今回私が作成したmmioptのver4.0の外部引数で入力できるパラメータにおいて、どの組み合わせを使うのが良さそうかを考察するため、私のPCで検証してみました。(詳細は実行環境の章を参照。)
思い切って、前回の記事の9倍の規模のTSPで検証してみましょう！
◆30000x30000.npz single no-2opt 固定（このパラメータはテキトーに決め打ちで実施。演算時間が早いので。）
| er※① | n※② | 距離 | 演算時間[s] |
|---|---|---|---|
| 0.00 | 29999 | 30351 | 4.7850 |
| 0.01 | 3137 | 30407 | 4.1949 |
| 0.05 | 125 | 30453 | 4.1029 |
| 0.10 | 31 | 30460 | 4.0309 |
| 0.30 | 3 | 30488 | 4.0099 |
| →やはり評価値の算出の許容誤差を増加させると、総距離が長いパスが出てきます。(^^; | |||
| 実用的には0.00～0.01を使うのが良さそうです |
※①er…許容誤差率です。
※②n…統計的なサンプルサイズnと同義。
nを求める計算式に、一様分布(1%～100%と決め打ち)の標準偏差を代入して作っています。
●一様分布の標準偏差 σ
https://bellcurve.jp/statistics/course/8013.html
1％～100％の一様分布の標準偏差は
σ=\sqrt{(1.00-0.01)^2/12}=0.2857884
(入力される行列一要素の値は相対値なので、
1%=1要素の持つ値域の幅/99 なので、値域により揺らぎますが、
ざっくりこう置けます。
1~100の自然数が入力されることを想定して作っています。)
●サンプルサイズ n
https://toukeigaku-jouhou.info/2018/01/23/how-to-calculate-samplesize/
n=(1.96*σ/許容できる誤差er)^2
n=int((1.96*0.2857884/er)^2)
プログラム上は、高速化のため、左端の[0]から、n個(対角要素は対象外)採取します。
実質、乱択アルゴリズムです。
◆30000x30000.npz er=0.01 固定
| 判断論理※ | 2-optあり/なし | 距離 | 演算時間[s] |
|---|---|---|---|
| single | no-2opt | 30407 | 4.1949 |
| single | 2opt | 30407 | 9.1305 |
| multi | no-2opt | 30367 | 221.0364 |
| multi | 2opt | 30367 | 228.9505 |
| ※解探索中に、複数の最大値/最小値に出くわしたときの判断論理 |
single… 元ノードの持つエッジが、複数の最小のエッジを持つ場合、強制的に若番ノードをその元ノードの解(先ノード)とする
multi…元ノードの持つエッジが、複数の最小のエッジを持つ場合、その中から、最大の評価値をもつノードを元ノードの解(先ノード)とする
no-2opt…解の2opt法の適用なし
2opt …解の2opt法の適用あり
multiだと時間のわりに、劇的に改善しないですね...
single/multiにおいて2-optを使用しても、結果として、変化なしでした。（そういうこともあるのか...というより、解の良さが最初のinitalizeとsingle/multiに依存する？）
ハイパーパラメータ考察結論
計算資源に余裕がない→single,er=0.01
計算資源に余裕がある→multi,er=0.00
no-2optで十分。
よって、実行コマンドは
計算資源に余裕がない場合は、
python msi_shortest4.0.py ファイル名.npz txt 0.01 single no2-opt
計算資源に余裕がある場合は、
python msi_shortest4.0.py ファイル名.npz txt 0.01 multi 2-opt
をおすすめします。(｀・ω・´)
社会インフラ（交通渋滞緩和や、ロボットなど）に組み込まれたら、すごく嬉しい...
実行環境
SW:
Python:Python 3.8.0
コマンドプロンプトより実行
HW:
OS:Windwos 10 Home
CPU:AMD Ryzen 7 2700 Eight-Core Processor 3.20GHz(OCなし,16論理コア)
DIMM:16.0GB 3200MHz
HDD:1TB
GPU:AMD Radeon R9 200
付録A (RaspberryPiに解かせてみた)
pi@raspberrypi:~/trunk/mmiopt $ python msinit_shortest1.0.py 10000x10000.npz txt 0.01 single no-2opt
start load file!
[[ 0 15 21 ... 7 14 75]
[56 0 94 ... 56 47 89]
[15 87 0 ... 68 88 56]
...
[ 5 10 79 ... 0 96 68]
[29 42 58 ... 20 0 56]
[16 76 49 ... 15 95 0]]
end load file!
('flags =', array([0, 0, 0, ..., 0, 0, 0]))
start max mean initalize.
er=0.01
n=3137
end max mean initialize.
start search min edges.
end search min edges.
start make ansedge from ansnode.
end make ansedge from ansnode.
('search end ansnode->', array([3596, 8, 197, ..., 9889, 9975, 9952]))
('Accumulation Distance = ', 10408.0)
('ansedge = ', array([ 1., 1., 1., ..., 15., 46., 47.]))
calculation elapsed time=8.08861899376
total elapsed time=13.6909029484
付録B (ソースコード)
# coding: UTF-8
# Maximum Mean Initialize and 2-opt.
# -> MMI-OPT
# this version is approximataly solve "shortest path" from TSP Input.
#
# This argolithm developed for solve TSP approximately fast.
# (TSP=Traveling Salesman Problem. )
#
# Confirmed operating environment.
# CentOS Linux release 7.2.1511 (Core)
# Python 2.7.5
# Python 3.8.0
# numpy 1.16.5
#
# How to use this software from command line?
# > python mmiopt_shortest4.0.py filename.npz txt 0.01 single no2-opt
# ^^^A solution is output as .npz
# ^^^^^^^^^^^^^int or float valuses square matrix
# of np.savez_compressed(filename.npz,array_=a)
# type
# Input:filename int or float
# Output:ansedge.npz or txt int
# ansnode.npz or txt int
# accumulation.npz or txt float
# developer "Ryota Tanaka"
import numpy as np
import csv
import sys
import time
import itertools as itr
def main():
argv = sys.argv
if len(argv) != 6:
print("File name may not be specified,or may contain extra options.")
exit()
print("start load file!")
data = np.load(argv[1])['array_']
print(data)
print("end load file!")
start = time.time()
searched = [data.shape[1]]
if data.shape[0] != data.shape[1]:
print("Please check if it is a square matrix.")
#array init
flags = np.zeros(data.shape[0],dtype='int32')
print("flags =",flags)
ansnode = np.array([],dtype='int32')
ansedge = np.array([],dtype='int32')
#ansedge = np.array([])
accumulation = 0 #累積距離の初期化
print("start max mean initalize.")
er=float(argv[3])
print("er="+str(er))
if er < 0.0 or 0.99 < er:
print("inputed error rate is out of range.")
exit()
elif er <=0.0:
n=data.shape[0]-1
else:
n=int((1.96*0.2857884/er)**2)
print("n="+str(n))
if n > data.shape[0]-1:
n=data.shape[0]-1
index,dataeva = calc_maxevaindex(data,n)
print("end max mean initialize.")
flags[index] = 100 #訪問済みにする
ansnode=np.append(ansnode,int(index)) #答えとなるノード番号をキューに入れる
num = data.shape[0] - 1
print("start search min edges.")
while num > 0:
beforeindex = index #ひとつ前のindexとして保持する
data[index]=np.add(data[index],flags)
if argv[4] == 'single':
index = np.argmin(data[index])
else:
#multi
xmin = np.min(data[index])
minindexes = [i for i,x in enumerate(data[index]) if x == xmin]
if len(minindexes) > 1:
maxdataevai = np.argmax(dataeva[minindexes])
index = minindexes[maxdataevai] #minindexesの内、評価値が最大のノード番号を取り出す
else:
index = minindexes[0]
flags[index] = 100 #訪問済みにする
ansnode=np.append(ansnode,int(index)) #答えとなるノード番号をキューに入れる
accumulation += data[beforeindex,index] #累積距離に加算する
num-=1
#print(num)
print("end search min edges.")
print("start make ansedge from ansnode.")
for i in range(len(ansnode)-1):
ansedge=np.append(ansedge,data[int(ansnode[i]),int(ansnode[i + 1])])
print("end make ansedge from ansnode.")
#ansedge = np.append(ansedge,data[int(ansnode[i]),int(ansnode[i + 1])])
#2-opt法による解のworst除去
if argv[5] == '2-opt':
print("start 2-opt")
num = data.shape[0] - 1
tempedge0=np.zeros(4)#従来のエッジ長集合
tempedge1=np.zeros(4)#swap後エッジ長集合
templist=np.zeros(6)
while num > 0:
worst_node_i = np.argmax(ansedge)#解の経路から、最長のエッジを出す
#次のループのために、現在見つかった最大値の要素を削除
ansedge = np.delete(ansedge,worst_node_i,axis=0)
worst_minnode=np.argmin(data[int(ansnode[worst_node_i])])
#◆交換なし
#交換元
worst_minnode_i_list=np.where(ansnode-worst_minnode==0)#array[]->intに変換は[0]を参照すればよい
worst_minnode_i=worst_minnode_i_list[0]
#print("worst_minnode_i");#print(worst_minnode_i)
worst_node_i_next = worst_node_i+1
#以下のif文は、worst_maxnode_iが解の先頭か、最後である場合、前がなかったり後がなかったりするので、その境界の処理
templist[1]=ansnode[worst_minnode_i]
if worst_minnode_i-1 > 0 :
templist[0]=ansnode[worst_minnode_i-1]
tempedge0[0]=data[int(templist[0]),int(templist[1])]
if worst_minnode_i+1 < len(ansnode)-1 :
templist[2]=ansnode[worst_minnode_i+1]
tempedge0[1]=data[int(templist[1]),int(templist[2])]
#交換先
templist[4]=ansnode[worst_node_i_next]
if worst_node_i-1 > 0 :
templist[3]=ansnode[worst_node_i_next-1]
tempedge0[2]=data[int(templist[3]),int(templist[4])]
if worst_node_i+1 < len(ansnode)-1 :
templist[5]=ansnode[worst_node_i_next+1]
tempedge0[3]=data[int(templist[4]),int(templist[5])]
#◆交換した場合
#交換元
templist[1]=ansnode[worst_node_i_next]#swaped
if worst_minnode_i-1 > 0 :
templist[0]=ansnode[worst_minnode_i-1]
tempedge1[0]=data[int(templist[0]),int(templist[1])]
if worst_minnode_i+1 < len(ansnode)-1 :
templist[2]=ansnode[worst_minnode_i+1]
tempedge1[1]=data[int(templist[1]),int(templist[2])]
#交換先
templist[4]=ansnode[worst_minnode_i]#swaped
if worst_node_i-1 > 0 :
templist[3]=ansnode[worst_node_i_next-1]
tempedge1[2]=data[int(templist[3]),int(templist[4])]
if worst_node_i+1 < len(ansnode)-1 :
templist[5]=ansnode[worst_node_i_next+1]
tempedge1[3]=data[int(templist[4]),int(templist[5])]
#print("sum(tempedge0)");#print(np.sum(tempedge0))
#print("sum(tempedge1)");#print(np.sum(tempedge1))
#距離判定
if(np.sum(tempedge1)<np.sum(tempedge0)):
#node swap
tempnode= ansnode[worst_node_i_next]
ansnode[worst_node_i_next]=ansnode[worst_minnode_i]
ansnode[worst_minnode_i]=tempnode
#decrement
num-=1
#while num end
print("end 2-opt")
print("search end ansnode->",ansnode)
ansedge =[]
for i in range(len(ansnode)-1):
ansedge = np.append(ansedge,data[int(ansnode[i]),int(ansnode[i + 1])])
accumulation=np.sum(ansedge)
print("Accumulation Distance = ",accumulation)
print("ansedge = ",ansedge)
if argv[2] == 'txt':
np.savetxt("ansnode.txt", ansnode,fmt ='%d')
np.savetxt("ansedge.txt", ansedge,fmt ='%d')
np.savetxt("accumulation.txt", [accumulation],fmt ='%f')
elif argv[2] == 'npz':
np.savez_compressed("ansnode.npz", array_=ansnode)
np.savez_compressed("ansedge.npz", array_=ansedge)
np.savez_compressed("accumulation.npz", array_=accumulation)
else:
print("!none output files!")
t = time.time() - start
print("calculation elapsed time="+str(t))
def calc_maxevaindex(data,n):
#dataeva = np.mean(data,axis=1)
maxj=data.shape[0]-1
#before half(row number 0 to n-1 -> sum(0 to n) (n+1 elements,include 0))
#print(data)
dataeva=[np.sum(data[0:n,:n+1],axis=1)]
#print("data[0:n,:n+1]="+str(data[0:n,:n+1]))
#after half(row number n to maxj -> sum(0 to n-1 (n elements))
#print("data[n:maxj+1,:n])"+str(data[n:maxj+1,:n]))
dataeva=np.append(dataeva,np.sum(data[n:maxj+1,:n],axis=1))
#print(dataeva)
#print("data[n:maxj+1,:n]="+str(data[n:maxj+1,:n])
#print("i=0 ="+str(dataeva))
#print(dataeva)
#i=1~maxj-1 (To get the line from which a diagonal element was removed.)
#print("maxj = "+str(maxj))
#print(str([np.mean(np.append(row[0:(i+1)],row[(i+1)+1:n+1],axis=0)) for i,row in enumerate(data[1:maxj])]))
#print(str(np.append(dataeva,[np.mean(np.append(row[0:(i+1)],row[(i+1)+1:n+1],axis=0)) for i,row in enumerate(data[1:maxj])],axis=0)))
#exit()
evamax_i=np.argmax(dataeva)
return evamax_i,dataeva
if __name__ == "__main__":
try:
start = time.time()
main() #main処理
t = time.time() - start
print("total elapsed time="+str(t))
#errorハンドリング
except IOError as error:
print("Make sure the files are on the same level")