ponder.ooo
🐍🐍🐍
laplacian growth
+193
-18
+6
-2
lib/log.py
+6
-2
lib/log.py
···
63
63
if file == "<string>" and source is not None:
64
64
file = inspect.getfile(source)
65
65
lines, first_line = inspect.getsourcelines(source)
66
-
line = lines[line_number-1].strip()
67
-
line_number = first_line + line_number - 1
66
+
if line_number < len(lines):
67
+
line = lines[line_number-1].strip()
68
+
line_number = first_line + line_number - 1
69
+
else:
70
+
# TODO fix this!!
71
+
_log(f"SNAKEPYT", "TRACING ERROR", mode="error", indent=self._indent+8)
68
72
file_end = file.split("/")[-1]
69
73
_log(f"in {file_end} {line_number}", line, mode="error", indent=self._indent+4)
70
74
return self
+11
main.py
+11
main.py
···
6
6
import sys
7
7
import os
8
8
import time
9
+
import shutil
9
10
from pathlib import Path
10
11
from argparse import ArgumentParser as ArgParser
11
12
from importlib import import_module, reload
···
19
20
#parser.add_argument("sketch", help="the sketch to run", type=str)
20
21
args = parser.parse_args()
21
22
23
+
snakepyt_version = (0, 0)
22
24
23
25
def establish_scheduler():
24
26
schedule = []
···
144
146
persistent_state = {}
145
147
persistent_hashes = {}
146
148
pyt_print("state flushed")
149
+
if command in ["exit", "quit", ":q", ",q"]:
150
+
pyt_print.blank().log("goodbye <3").blank()
151
+
repl_continue = False
152
+
continue
147
153
if command == "run":
148
154
sketch_name, remainder = lsnap(remainder)
149
155
···
182
188
run_dir = time.strftime(f"%d.%m.%Y/{sketch_name}_t%H.%M.%S")
183
189
sketch.__dict__["run_dir"] = run_dir
184
190
Path("out/" + run_dir).mkdir(parents=True, exist_ok=True)
191
+
192
+
shutil.copy(sketch.__file__, f"out/{run_dir}/{sketch_name}.py")
193
+
194
+
with open(f"out/{run_dir}/.snakepyt", "w") as metadata:
195
+
metadata.write(f"snakepyt version {snakepyt_version[0]}.{snakepyt_version[1]}\n")
185
196
186
197
if hasattr(sketch, "init"):
187
198
if hasattr(sketch, "final"):
+142
sketch/laplace.py
+142
sketch/laplace.py
···
1
+
# laplacian growth
2
+
# see "Fast Simulation of Laplacian Growth" by Theodore Kim, Jason Sewall, Avneesh Sud and Ming C. Lin
3
+
4
+
import itertools
5
+
from math import tau, log
6
+
7
+
import torch
8
+
9
+
from lib.util import load_image_tensor, save, msave
10
+
11
+
12
+
13
+
def init():
14
+
torch.set_default_device("cuda")
15
+
16
+
dims = 2
17
+
#scale = 2 ** 10
18
+
#h = scale
19
+
#w = scale
20
+
margin = 2**4
21
+
22
+
physics_scale = 0.5
23
+
eta = 10.7 # fractal dimension or something
24
+
25
+
max_points = 2 ** 13
26
+
27
+
t_real = torch.half
28
+
29
+
point_positions = torch.zeros([max_points, dims], dtype=t_real)
30
+
point_values = torch.zeros([max_points], dtype=t_real)
31
+
is_charge = torch.zeros([max_points], dtype=torch.bool)
32
+
do_neighbors = torch.zeros([max_points], dtype=torch.bool)
33
+
is_candidate = torch.zeros_like(is_charge)
34
+
is_free = torch.ones_like(is_charge)
35
+
36
+
neighbors = torch.tensor([v for v in itertools.product([-1,0,1], repeat=dims) if any(v)], dtype=t_real)
37
+
38
+
inner_radius = physics_scale / 2
39
+
40
+
def insert_charge(position, value, was_candidate=True, no_neighbors=False):
41
+
if was_candidate: # index could probably be provided by the caller for efficiency
42
+
match_coords = point_positions == position.expand(point_positions.shape)
43
+
match_pos = match_coords.prod(dim=1)
44
+
index = torch.nonzero(is_candidate & match_pos)[0]
45
+
is_candidate[index] = False
46
+
is_charge[index] = True
47
+
do_neighbors[index] = not no_neighbors
48
+
point_values[index] = value
49
+
else:
50
+
# create the charge
51
+
free_index = torch.nonzero(is_free)[0] # maybe not ideal
52
+
is_charge[free_index] = True
53
+
is_free[free_index] = False
54
+
do_neighbors[free_index] = not no_neighbors
55
+
point_positions[free_index] = position
56
+
point_values[free_index] = value
57
+
58
+
# update existing candidate potentials (eqn 11)
59
+
candidate_pos = point_positions[is_candidate]
60
+
dist = (candidate_pos - position.expand(candidate_pos.shape)).norm(p=2, dim=1)
61
+
point_values[is_candidate] += 1 - inner_radius/dist
62
+
63
+
if no_neighbors:
64
+
return
65
+
66
+
# add new candidates
67
+
charge_pos = point_positions[is_charge & do_neighbors]
68
+
for charge_index in range(charge_pos.shape[0]):
69
+
neighborhood = neighbors + charge_pos[charge_index].expand(neighbors.shape)
70
+
for neighbor_index in range(neighborhood.shape[0]):
71
+
maybe_insert_candidate(neighborhood[neighbor_index])
72
+
73
+
def maybe_insert_candidate(position):
74
+
if torch.any((point_positions == position.expand(point_positions.shape)).prod(dim=1)):
75
+
return
76
+
77
+
free_index = torch.nonzero(is_free)[0] # maybe not ideal
78
+
is_candidate[free_index] = True
79
+
is_free[free_index] = False
80
+
point_positions[free_index] = position
81
+
charge_pos = point_positions[is_charge]
82
+
charge_val = point_values[is_charge]
83
+
dist = (charge_pos - position.expand(charge_pos.shape)).norm(p=2, dim=1)
84
+
point_values[free_index] = (1 - inner_radius/dist).sum() # eqn 10
85
+
86
+
def select_candidate(scale):
87
+
try:
88
+
candidate_val = point_values[is_candidate]
89
+
candidate_pos = point_positions[is_candidate]
90
+
max_val = candidate_val.max()
91
+
min_val = candidate_val.min()
92
+
normalized_val = (candidate_val - min_val) / (max_val - min_val) # eqn 13
93
+
val_pow = normalized_val ** eta
94
+
selection_prob = val_pow / val_pow.sum()
95
+
96
+
prob_sum = selection_prob.cumsum(dim=0)
97
+
98
+
choices = []
99
+
for _ in range(scale):
100
+
r = torch.rand(1, dtype=t_real)
101
+
choice = torch.nonzero(prob_sum > r)[0]
102
+
if choice not in choices:
103
+
insert_charge(candidate_pos[choice], 1)
104
+
choices.append(choice)
105
+
except KeyboardInterrupt:
106
+
raise
107
+
except:
108
+
print(f"selection failed")
109
+
110
+
111
+
schedule(run, None)
112
+
113
+
114
+
def run():
115
+
116
+
# set initial conditions
117
+
insert_charge(torch.tensor([0]*dims, dtype=t_real), 1, was_candidate=False)
118
+
for n in range(-40,40):
119
+
insert_charge(torch.tensor([2, n], dtype=t_real), 1, was_candidate=False, no_neighbors=True)
120
+
121
+
schedule(main_loop, None)
122
+
123
+
def main_loop():
124
+
for iteration in range(10000):
125
+
select_candidate(int(log(iteration * 100 + 10)))
126
+
if iteration % 10 == 0:
127
+
final_charges = point_positions[is_charge].clone()
128
+
129
+
for d in range(dims):
130
+
final_charges[:,d] -= final_charges[:,d].min()
131
+
final_charges[:,d] += margin
132
+
133
+
h = int(final_charges[:,0].max().item() + margin)
134
+
w = int(final_charges[:,1].max().item() + margin)
135
+
136
+
indices = final_charges.long()
137
+
138
+
canvas = torch.ones([h,w], dtype=t_real)
139
+
canvas[(indices[:,0], indices[:,1])] = 0
140
+
141
+
msave(canvas, f"{run_dir}/{iteration:06d}")
142
+
+34
-16
sketch/plume/plume_c.py
+34
-16
sketch/plume/plume_c.py
···
22
22
random_range = tau #/ 50
23
23
24
24
# 2 ** 9 = 512; 10 => 1024; 11 => 2048
25
-
scale_power = 12
25
+
scale_power = 13
26
26
scale = 2 ** scale_power
27
27
28
28
···
39
39
out[:,dim] *= scale
40
40
return out
41
41
42
+
def bad_proj(u, v):
43
+
dot = (u * v).sum(dim=1)
44
+
scale = (dot / v.norm(p=2, dim=1))
45
+
out = v.clone()
46
+
for dim in range(out.shape[1]):
47
+
out[:,dim] *= scale
48
+
return out
49
+
42
50
def proj_shift(a, b):
43
51
return a + proj(b, a)
52
+
53
+
def bad_proj_shift(a, b):
54
+
return a + bad_proj(b, a)
44
55
45
56
def get_transformation(direction, flow, ebb, rescaling):
46
57
def transformation(p):
58
+
d = direction
59
+
d[:,0] += torch.rand_like(direction[:,0]) * 0.0001
47
60
if flow == 0:
48
61
result = p - direction * ebb
49
62
else:
50
-
result = proj_shift(p, direction * flow) - direction * ebb
63
+
result = proj_shift(p, d * flow)
64
+
result -= direction * ebb #proj_shift(result, direction * ebb)
51
65
if rescaling:
52
66
result /= result.norm(p=2,dim=1).max()
53
67
return result
···
61
75
return math.pow(b/a, t)*a
62
76
63
77
def per_t(t):
64
-
origin = 0, -1
78
+
origin = 0.0625, -0.15
65
79
66
80
#s = eerp(1000, 0.001, t)
67
-
s = 2
68
-
x_range = 1.2 * s
69
-
y_range = 1.2 * s
81
+
s = 0.25
82
+
x_range = 0.5 * s
83
+
y_range = 1 * s
70
84
71
85
span = x_range, y_range
72
86
mapping = map_space(origin, span, zooms, stretch, scale)
···
82
96
direction = torch.ones_like(p_positions)
83
97
direction[:,0] = 0
84
98
85
-
iterations = 100
99
+
iterations = 1301
86
100
#show_frame = lambda i: i == iterations - 1
87
-
show_frame = lambda i: i % 10 == 0#True
101
+
show_frame = lambda i: i % 100 == 0#True
88
102
#show_frame = lambda i: True
89
103
90
-
ebb = 0.804#0.9000469530469531 - 0.001 + 0.002 * t
104
+
#ebb = 0.804#0.9000469530469531 - 0.001 + 0.002 * t
105
+
#ebb = 1 + 0.4 * t
106
+
#ebb = 0.4
107
+
ebb = 0.83314# + (0.001 * t)
91
108
#ebb = (0.5 + 0.5 * t)#0.804#482 / 600
92
-
#ebb =0.9000469530469531 #0.9000595404595405
109
+
#ebb = 0.0 + 0.01 * t #0.9000469530469531 #+ 0.00001 * t #0.9000595404595405
93
110
flow = 1
94
111
rescaling = False
95
112
···
122
139
#diff /= diff.max()
123
140
124
141
if show_frame(iteration):
142
+
#diff_avg = diff
125
143
frame_index[0] += 1
126
144
print(f"{frame_index[0]}: ebb={ebb}")
127
145
128
146
#pos_reshape = p_positions.reshape((h,w,2))
129
-
#scratch[:,:,0] = pos_reshape[:,:,0]
130
-
#scratch[:,:,2] = pos_reshape[:,:,1]
147
+
#scratch[:,:,0] = pos_reshape[:,:,0] + 0.5
148
+
#scratch[:,:,1] = pos_reshape[:,:,1] + 0.5
131
149
132
150
#scratch += 0.5
133
151
#scratch *= 2
134
-
scratch[:,:,0] = diff_avg.nan_to_num().reshape((h,w))
152
+
scratch[:,:,2] = diff_avg.nan_to_num().reshape((h,w))
135
153
#scratch[:,:,2] = diff_avg.nan_to_num().reshape((h,w))
136
154
137
155
scratch[scratch < 0] = 0
···
139
157
#scratch[:,:,0] = scratch[:,:,0].pow(2)
140
158
#scratch[:,:,2] = scratch[:,:,2].pow(0.5)
141
159
142
-
scratch[:,:,0] -= scratch[:,:,0].mean()
143
-
scratch[:,:,0] /= scratch[:,:,0].std() * 6
160
+
scratch[:,:,2] -= scratch[:,:,2].mean()
161
+
scratch[:,:,2] /= scratch[:,:,2].std() * 6
144
162
#scratch[:,:,2] /= scratch[:,:,2].max()
145
163
146
-
scratch[:,:,0] += 0.5
164
+
scratch[:,:,2] += 0.5
147
165
148
166
#scratch[:,:,1] = diff.reshape((h,w)) / diff.max()
149
167