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

```txt:動作環境
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](http://epubs.siam.org/doi/abs/10.1137/140960463?journalCode=sjnaam)
参考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値取得関数のチェック](https://qiita.com/7of9/items/c94a4b9fd85a972c00fc)
においてtriangular van der Corput列に用いるtriangleへの分岐経路を取得することができるようになった。

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


### code

```py: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

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/8ab1b2ed-b111-55b3-00e1-7374162d8ed4.png)

#### NUM_SEQ_NO = 4**2

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/cfd1785e-db7c-1de6-891f-4670d409d890.png)

#### NUM_SEQ_NO = 4**3

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/41f1745e-124c-d14a-62c8-1f2cc0882454.png)

#### NUM_SEQ_NO = 4**4

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/e803be4f-9fe2-d711-058a-975d61b6acd4.png)

#### NUM_SEQ_NO = 4**5

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/7108e607-f982-a6a6-9387-a2d09fd4530a.png)

### 備考

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

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

![qiita.png](https://qiita-image-store.s3.amazonaws.com/0/32870/67d84bb4-62a9-7777-2bee-d084c04714ad.png)

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