QMC (Quasi-Monte Carlo Methods)で使うHalton SequenceのPython実装
v0.2
http://qiita.com/7of9/items/5e61943bcb898521e1c6
のC実装をPython用にしただけ。
def Halton_sequence(i0):
xbase = 2
ybase = 3
invxbase = 1.0 / xbase
facx = 1.0 / xbase
invybase = 1.0 / ybase
facy = 1.0 / ybase
inp = i0
x0 = 0.0
while inp > 0:
x0 = x0 + (inp % xbase) * invxbase
inp = inp / xbase
invxbase = invxbase * facx
inp = i0
y0 = 0.0
while inp > 0:
y0 = y0 + (inp % ybase) * invybase
inp = inp / ybase
invybase = invybase * facy
return x0, y0
for i0 in range(10):
x0,y0 = Halton_sequence(i0)
print x0,y0
結果
0.0 0.0
0.5 0.333333333333
0.25 0.666666666667
0.75 0.111111111111
0.125 0.444444444444
0.625 0.777777777778
0.375 0.222222222222
0.875 0.555555555556
0.0625 0.888888888889
0.5625 0.037037037037
見覚えのある数値が出ているので多分あっているかと。
時間ができたときに整理するかも。
v0.3
(追記 2017/03/20)
importできるように変更した。
static methodにしている。
参考 http://stackoverflow.com/questions/735975/static-methods-in-python
UtilHaltonSequence.py
'''
v0.3 Mar. 20, 2017
- change to a static function in a Class
v0.2 Oct. 22, 2016
- implemented in Python
v0.1 Mar., 2005 or so
- implemented in C
'''
# codingrule:PEP8
class CHaltonSequence:
@staticmethod
def calc_Halton_sequence(index):
XBASE = 2
YBASE = 3
inv_xbase = 1.0 / XBASE
fac_x = 1.0 / XBASE
inv_ybase = 1.0 / YBASE
fac_y = 1.0 / YBASE
inp = index
xwrk = 0.0
while inp > 0:
xwrk = xwrk + (inp % XBASE) * inv_xbase
inp = inp / XBASE
inv_xbase = inv_xbase * fac_x
inp = index
ywrk = 0.0
while inp > 0:
ywrk = ywrk + (inp % YBASE) * inv_ybase
inp = inp / YBASE
inv_ybase = inv_ybase * fac_y
return xwrk, ywrk
'''
Usage example:
from UtilHaltonSequence import CHaltonSequence
for idx in range(0, 11):
x0, y0 = CHaltonSequence.calc_Halton_sequence(index=idx)
print(x0, y0)
'''
test_haltonSeq.py
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from UtilHaltonSequence import CHaltonSequence
for idx in range(0, 11):
x0, y0 = CHaltonSequence.calc_Halton_sequence(index=idx)
print(x0, y0)
実行結果
$ python test_haltonSeq.py
0.0 0.0
0.5 0.333333333333
0.25 0.666666666667
0.75 0.111111111111
0.125 0.444444444444
0.625 0.777777777778
0.375 0.222222222222
0.875 0.555555555556
0.0625 0.888888888889
0.5625 0.037037037037
0.3125 0.37037037037