这篇文章发表于 1162 天前，可能其部分内容已经发生变化，如有疑问可询问作者。

突然感觉很需要去学一学这东西，以后可能自己写地理论文的时候会用到，于是我就入坑了……

taichi编个程序门槛挺高的，要一定的数理基础——然而我快忘完了，其次还有点不适应。结果它的Hello World程序看了半天，心态小崩。

分形图像

Julia Set

如果你兴冲冲地去看官方提供的Hello World程序的话，很可能会被泼冷水——主要困扰我的是其主循环中的i变量。后来发现其实这个变量是控制图形能够不断变化的。这里建议先看这篇文章，他用Pascal写的程序，他的这个图就是固定的（相当于i变量不变的情况）。

对照他的程序，我写了taichi版本的程序与其功能基本一致。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import taichi as ti

ti.init(arch=ti.gpu)

pixels = ti.field(dtype=float, shape=(600, 400))
c = ti.Vector([-0.75, 0])


@ti.func
def complex_sqr(z): # 复数平方
    return ti.Vector([z[0]**2 - z[1]**2, 2 * z[0] * z[1]])


@ti.kernel
def paint():
    for i, j in pixels:
        z = ti.Vector([i / 200, j / 200])
        iterations = 0
        while z.norm() < 20 and iterations < 200:
            z = complex_sqr(z) + c
            iterations += 1
        pixels[i, j] = 1 - iterations * 0.01


gui = ti.GUI("Julia Set", res=(600, 400))
while True:
    paint()
    gui.set_image(pixels)
    gui.show()

相信运行过这个程序的你肯定能很快明白那个Hello World程序的意思了。

juliaset

Sierpinski Triangle

接下来我们尝试写其他的一些分形图形来熟悉操作，这里选择了Sierpinski三角形。

对照这篇文章中的Pascal程序（他使用了Chaos Game的随机构造方法），taichi版本的如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import taichi as ti

ti.init(arch=ti.gpu)

width, height = 600, 400
pixels = ti.field(dtype=float, shape=(width, height))

N = 3
x = ti.Vector.field(2, float, N)
random_x = ti.field(float, ())
random_y = ti.field(float, ())


@ti.kernel
def init():
    """
    for i in x:
        x[i] = ti.Vector([ti.random(),
                          ti.random()])  # ti.random() belongs [0,1], and it does not have random.seed
    """
    x[0] = ti.Vector([0.9, 0.9])
    x[1] = ti.Vector([0.2, 0.4])
    x[2] = ti.Vector([0.5, 0.2])
    pixels.fill(1)


init()
random_x[None], random_y[None] = x[0][0], x[0][1]


@ti.func
def next_point(rand_x, rand_y):
    r = ti.random() * 3
    new_x = (x[r * 10 // 10][0] + rand_x) / 2
    new_y = (x[r * 10 // 10][1] + rand_y) / 2
    pixels[int(new_x * width), int(new_y * height)] = 0
    return ti.Vector([new_x, new_y])


@ti.kernel
def paint():
    tmp = next_point(random_x[None], random_y[None])
    random_x[None], random_y[None] = tmp[0], tmp[1]


gui = ti.GUI("Sierpinski Triangle", res=(width, height))
while True:
    paint()
    gui.set_image(pixels)
    gui.show()

效果就是这样子啦。

sierpinski

Barnsley Fern

看到论坛上有人写了一个BarnsleyFern巴恩斯利蕨的动图，于是自己也编一个。

这篇文章可以用来参考，但我参考的是wikipedia上的内容，对taichi而言，其代码完全可以变得更优雅。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import taichi as ti
ti.init(arch=ti.gpu)

width, height = 600, 400
pixels = ti.field(dtype=float, shape=(width, height))

x = ti.field(float, ())
y = ti.field(float, ())
ret = ti.Vector.field(2, float, shape=())
x[None], y[None] = 0, 0

A1 = ti.Matrix([[0, 0], [0, 0.16]])
A2 = ti.Matrix([[0.85, 0.04], [-0.04, 0.85]])
A3 = ti.Matrix([[0.2, -0.26], [0.23, 0.22]])
A4 = ti.Matrix([[-0.15, 0.28], [0.26, 0.24]])
C1 = ti.Vector([0, 0])
C2 = ti.Vector([0, 1.6])
C3 = ti.Vector([0, 1.6])
C4 = ti.Vector([0, 0.44])


@ti.func
def next_point(in_x, in_y):
    r = ti.random()
    new_x, new_y = int((in_x / 6 + 0.5) * width), int(in_y * height / 10.5)
    if r < 0.01:
        pixels[new_x, new_y] = 0.1
        ret[None] = A1 @ ti.Vector([in_x, in_y]) + C1
    elif r < 0.86:
        pixels[new_x, new_y] = 0.4
        ret[None] = A2 @ ti.Vector([in_x, in_y]) + C2
    elif r < 0.93:
        pixels[new_x, new_y] = 0.6
        ret[None] = A3 @ ti.Vector([in_x, in_y]) + C3
    else:
        pixels[new_x, new_y] = 0.9
        ret[None] = A4 @ ti.Vector([in_x, in_y]) + C4
    return ret[None]


@ti.kernel
def paint():
    tmp = next_point(x[None], y[None])
    x[None], y[None] = tmp[0], tmp[1]


gui = ti.GUI("Barnsley Fern", res=(width, height))

while True:
    paint()
    gui.set_image(pixels)
    gui.show()

于是就产生了这样的静态图片。

barnsleyfern

如果您理解了那个Hello World程序，那么想要让它动起来很简单。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
import taichi as ti
import math

ti.init(arch=ti.gpu)

width, height = 600, 400
pixels = ti.field(dtype=float, shape=(width, height))

x = ti.field(float, ())
y = ti.field(float, ())
ret = ti.Vector.field(2, float, shape=())

A1 = ti.Matrix([[0, 0], [0, 0.16]])
A2 = ti.Matrix([[0.85, 0.04], [-0.04, 0.85]])
A3 = ti.Matrix([[0.2, -0.26], [0.23, 0.22]])
A4 = ti.Matrix([[-0.15, 0.28], [0.26, 0.24]])
C1 = ti.Vector([0, 0])
C2 = ti.Vector([0, 1.6])
C3 = ti.Vector([0, 1.6])
C4 = ti.Vector([0, 0.44])


@ti.func
def next_point(in_x, in_y, t):
    r = ti.random()
    new_x, new_y = int((in_x / 6 + 0.5) * width), int(in_y * height / 10.5)
    if r < 0.01:
        pixels[new_x, new_y] = 0.1
        ret[None] = A1 @ ti.Vector([in_x, in_y]) + C1
    elif r < 0.86:
        pixels[new_x, new_y] = 0.4
        a = ti.Matrix([[0.85, 0.04 + t], [-0.04 - t, 0.85]])
        ret[None] = a @ ti.Vector([in_x, in_y]) + C2
    elif r < 0.93:
        pixels[new_x, new_y] = 0.6
        ret[None] = A3 @ ti.Vector([in_x, in_y]) + C3
    else:
        pixels[new_x, new_y] = 0.9
        ret[None] = A4 @ ti.Vector([in_x, in_y]) + C4
    return ret[None]


@ti.kernel
def paint(t: float):
    for i in range(0, 100000):
        x[None] = 0
        y[None] = 0
        for j in range(0, 40):
            tmp = next_point(x[None], y[None], t)
            x[None], y[None] = tmp[0], tmp[1]


gui = ti.GUI("Barnsley Fern", res=(width, height))
timer = 0

while True:
    paint(math.sin(timer) * 0.03)
    timer += 0.05
    gui.set_image(pixels)
    gui.show()
    pixels.fill(0)

barnsleyfern

引力

其实本来是想利用taichi写一个地理上引力模型的动图的，后来发现自己原来一直对这个概念有着错误的理解，那就在这里把这个迷思给破解了吧。

原本以为地理上的引力模型和物理上的引力模型基本上是同一个东西，但直到今天我才意识到两者虽然公式尽管一模一样却在意义上有着天壤之别：举个例子，平面中有两个固定的重量相同的圆球A和B，然后空间中又冒出一个初动能为零的重量相同的小球，试问小球最终会落到哪里？（从地理角度提的话，就是研究区域有两大商场，问该空间中的消费者会选择哪一个。）

显然最终落到哪里肯定和小球冒出的位置有关。但接下来地理和物理的思考视角就不同了。

地理上的引力模型说确切一点其实涉及到的是较为宏观的指标（像吸引到的零售额或是小球被吸引的概率，前者是赖利提出的，后者是赫夫的概率引力模型）的比值。另外，个人认为吸引力的大小的比较（谁吸引力大，小球就去谁那里）其实也可以作为一个简单粗暴的判定标准（之后的运动过程或是轨迹则不关心）。

但物理学肯定会认为谁吸引力大就去谁那里这种想法是错误的，因为这不是一个力的简单比较。小球受到吸引后会运动，在运动过程中圆球A和B的吸引力也会变化，不可能像地理上的引力模型那样靠算个力就能直接得出结果。

地理

地理上的没啥可多说的。像赖利的理论压根不是从单个消费者的角度考虑的；而赫夫的理论涉及到了概率，其在空间上的分布不会是单一的简洁的形态。但如果仅考虑吸引力大小，那么我们就能利用“二力平衡”找到边界，最后边界形成一个（阿波罗尼斯）圆（当然，如果两者引力相当那就是一条直线的特殊情况）。

物理

个人感觉物理上的这个问题似乎并不简单，至少还和A，B和小球的半径有关，那就假设A和B半径一样吧，而小球的半径就不计了。另外，当小球离A和B距离足够小的时候判定为“落”在A和B球上。

由于自己目前物理能力感人，我只能借助CPU了。这部分的代码我修改自taichi的example中的ad_gravity.py。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import taichi as ti
ti.init()

dt, N = 1e-5, 400

x = ti.Vector.field(2, float, N, needs_grad=True)  # position of particles
status = ti.field(int, N)  # if hit A, status[i] = 1; hit B, 2; else 0
origin_x = ti.Vector.field(2, float, N)  # origin position of particles
fixed = ti.Vector.field(2, float, 2)  # position of ball A and B

v = ti.Vector.field(2, float, N)  # velocity of particles
U = ti.field(float, (), needs_grad=True)  # potential energy

hit_A = ti.Vector.field(2, float, N)
hit_B = ti.Vector.field(2, float, N)
num_hit_A = ti.field(int, ())
num_hit_B = ti.field(int, ())
hit_A.fill(0)
hit_B.fill(0)
num_hit_A[None], num_hit_B[None] = 0, 0
status.fill(0)


@ti.kernel
def compute_U():
    for i in range(N):
        for j in range(2):
            r = x[i] - fixed[j]
            U[None] += -1 / r.norm(1e-3)  # U += -1 / |r|


@ti.kernel
def advance():
    for i in x:
        if status[i] == 0:
            v[i] += dt * -x.grad[i]  # dv/dt = -dU/dx
            x[i] += dt * v[i]  # dx/dt = v
        else:
            v[i] = ti.Vector([0, 0])


@ti.kernel
def judge():
    for i in x:
        if status[i] == 0:
            if (x[i] - fixed[0]).norm() < 0.005:
                hit_A[num_hit_A[None]] = origin_x[i]
                num_hit_A[None] += 1
                status[i] = 1
            elif (x[i] - fixed[1]).norm() < 0.005:
                hit_B[num_hit_B[None]] = origin_x[i]
                num_hit_B[None] += 1
                status[i] = 2


def substep():
    with ti.Tape(U):
        compute_U()
    advance()
    judge()


@ti.kernel
def init():
    fixed[0], fixed[1] = [0.35, 0.5], [0.65, 0.5]
    for i in x:
        x[i] = [ti.random(), ti.random()]
        origin_x[i] = x[i]


init()
gui = ti.GUI('Solution')
while gui.running:
    for i in range(50):
        substep()
    #print('U = ', U[None])
    gui.circles(origin_x.to_numpy(), radius=3, color=0x869a88)
    gui.circles(fixed.to_numpy(), radius=5, color=0xDDDDDD)
    gui.circles(hit_B.to_numpy(), radius=4, color=0xe4f35b)
    gui.circles(hit_A.to_numpy(), radius=4, color=0xf90c4b)
    gui.circles(x.to_numpy(), radius=3)
    gui.show()

当N=400的时候以上代码运行效果大致如下（越看越像东方）：

gravitation

上面两个对称的大白点就是固定球体A和B，而周围其他生成的随机点一旦过于接近球体A和B就认为是“落”在其上，并对其原始位置通过颜色标记来完成区分（红色是“落”在了左边A球上的小球的原始位置，而黄色是“落”在了右边B球上的小球的原始位置，灰色是尚在运动中的小球的原始位置）——这样我们只需要增加小球的数量就能够发现一些有趣的事情。

当然这么点小球肯定是不够的，于是我们令N=8000，运行一段时间后得到了如下的结果：

gravity

从这张让人摸不着头脑的图片来看，即便一些情况被大大简化了，这个物理题依然不简单。