geometry > sub-triangles > Triangular van der Corput列 > v0.2 > Matplotlibでの描画をしてみた

15699 ワード
matplotlib borgWarp geometry QMC matplotlib テキストリンク
動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)

参考1: Basu and Owen, SIAM J. NUMER. ANAL. 53, 743-761, 2015
参考2: http://statweb.stanford.edu/~owen/pubtalks/mcm2017-annotated.pdf
(p15あたりから)

QMC(Quasi Monte Carlo)法で使ったことがあるvan der Corput列。
Basu and Owen(2015)によって、Triangular van der Corput列が研究されていることを知った。

手順の概要は参考1に記載されている。

インデックスに基づき三角形を４当分していき、該当の三角形のcentroidを結果として使う。

関連

v0.1: https://qiita.com/7of9/items/ed5d8146aa03a5be2221

v0.2

geometry > 連番から分岐経路番号を取得する実装 (上部:4分岐, 最下部のみ:3分岐)> v0.5: 直値でのd値取得関数のチェック
においてtriangular van der Corput列に用いるtriangleへの分岐経路を取得することができるようになった。

それを元にvan der Corput列を取得してみた。
その結果を図示してみる。

code

triangularVanDerCorput_171017.ipynb
%matplotlib inline

# v0.2 Oct. 21, 2017
#  - increase figure size
#  - change loop structure to use values from convert_to_dvalues()
#  - merge [subtriangle_171014.py: v0.2]
# v0.1 Oct. 17, 2017
#  - branched from [test_scatter_170812.ipynb]

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

def get_centroid(vertices):
    inA, inB, inC = vertices
    return (inA + inB + inC) / 3


def get_subvertices(outer, didx):
    # outer: vertices of original triangle
    # didx: [0..3]
    #
    # Reference: Basu and Owen, SIAM J. NUMER. ANAL. 53, 743-761, 2015
    #
    inA, inB, inC = outer
    if didx >= 4:
        return 0, 0, 0  # error
    if didx == 0:
        resA = (inB + inC) / 2
        resB = (inA + inC) / 2
        resC = (inA + inB) / 2
    if didx == 1:
        resA = inA
        resB = (inA + inB) / 2
        resC = (inA + inC) / 2
    if didx == 2:
        resA = (inB + inA) / 2
        resB = inB
        resC = (inB + inC) / 2
    if didx == 3:
        resA = (inC + inA) / 2
        resB = (inC + inB) / 2
        resC = inC
    return resA, resB, resC


def convert_to_dvalues(decimal, base):
    if decimal < base:
        return [decimal]
    res = [(decimal - 1) % (base-1) + 1]
    while decimal >= base:
        res += [(((decimal - base) * base // (base-1)) // base) % base]
        decimal //= base
    res.reverse()
    return res


# regular triangle
VTXA = np.array([1, 0])
VTXB = np.array([-1, 0])
VTXC = np.array([0, 2])

MAX_DEPTH = 5
DVAL_BASE = 4

NUM_SEQ_NO = 4**5   # should be 4**N

xs = []
ys = []
for seqno in range(NUM_SEQ_NO):
    dvals = convert_to_dvalues(seqno, DVAL_BASE)

    resvtx = [VTXA, VTXB, VTXC]
    for dpos in range(len(dvals) - 1):
        resvtx = get_subvertices(resvtx, didx=dvals[dpos])
    resvtx = get_subvertices(resvtx, didx=dvals[-1])
    cntrd = get_centroid(resvtx)
    xs += [cntrd[0]]
    ys += [cntrd[1]]

color = [10 for idx in range(len(xs))]

# 2. draw data
size = 5  # arbitrary
fig = plt.figure(dpi=110)
plt.scatter(xs, ys, size, color, cmap=cm.binary)
plt.colorbar()

結果

NUM_SEQ_NO = 4**1

NUM_SEQ_NO = 4**2

NUM_SEQ_NO = 4**3

NUM_SEQ_NO = 4**4

NUM_SEQ_NO = 4**5

備考

「NUM_SEQ_NO=4の乗数」の場合は参考2の結果と同じものが得られたと思う。

一方で「NUM_SEQ_NO=32」の結果は以下のようになり、参考２とは異なる。

「4の乗数でない値」までのtriangular van der Corput列を用いた積分の精度がどうなるのか（Low Discrepancyはどうなるか）については未消化。