LoginSignup
1

More than 5 years have passed since last update.

Python > QMC(Halton_Sequence)の実装 > v0.2, v0.3

Last updated at Posted at 2016-10-12

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1