geometry > sub-triangles > Triangular van der Corput列 > v0.1 > for文のネスト

動作環境

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を結果として使う。

code v0.1

ネスト方式での実装

subtriangle_171014.py

import numpy as np

# on Python 3.5.2

# v0.1 Oct. 14, 2017
#   - add get_centroid()
#   - loop until two depths
#   - add get_subvertices()


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

# regular triangle
pA = np.array([1, 0])
pB = np.array([-1, 0])
pC = np.array([0, 2])

# loop
for d1 in range(4):
    s1A, s1B, s1C = get_subvertices([pA, pB, pC], didx=d1)
    cntrd = get_centroid([s1A, s1B, s1C])
    print(cntrd)
    for d2 in range(4):
        if d2 == 0:
            continue
        s2A, s2B, s2C = get_subvertices([s1A, s1B, s1C], didx=d2)
        cntrd = get_centroid([s2A, s2B, s2C])
        print(cntrd)

run

$ python3 subtriangle_171014.py 
[ 0.          0.66666667]
[-0.25        0.83333333]
[ 0.25        0.83333333]
[ 0.          0.33333333]
[ 0.5         0.33333333]
[ 0.75        0.16666667]
[ 0.25        0.16666667]
[ 0.5         0.66666667]
[-0.5         0.33333333]
[-0.25        0.16666667]
[-0.75        0.16666667]
[-0.5         0.66666667]
[ 0.          1.33333333]
[ 0.25        1.16666667]
[-0.25        1.16666667]
[ 0.          1.66666667]

本来の式

本来は以下の式に基づく。

T(d1, d2) = (T(d1))(d2)

d2は0以外。

与えられたインデックスを元にd1, d2の値を求めて、逐次T()を計算していけばいいはず。