1	# bend-visgraph's library: ts2vg's options and a series in, the graph out,
2	# with no def that relies on @unsafe. Its graphs are ts2vg 1.2.4's, in F32.
3	# parallel.bend builds the same graphs on every thread, and relies on Base's
4	# @unsafe Array.fork; a program that does not import it passes
5	# `bend --check-only`.
6	#
7	#   import bend-visgraph@2.0.0.0/visgraph.bend as VG
8	#
9	# Its interface is the types at the top, inf and neg_inf, and the defs under
10	# "The library" at the bottom; every other def serves them.
11	import Base
12	
13	type Visibility is Data:
14	  Natural{}
15	  Horizontal{}
16	
17	type Direction is Data:
18	  Undirected{}
19	  LeftToRight{}
20	  TopToBottom{}
21	
22	type Weight is Data:
23	  Unweighted{}
24	  Distance{}
25	  SqDistance{}
26	  VDistance{}
27	  AbsVDistance{}
28	  HDistance{}
29	  AbsHDistance{}
30	  Slope{}
31	  AbsSlope{}
32	  Angle{}
33	  AbsAngle{}
34	  NumPenetrations{}
35	
36	# The x coordinate of each value: its index, or the value at the same
37	# index of an array.
38	type Xs is Type:
39	  Index{}
40	  At{xs: Array<F32>}
41	
42	# The edges, as a rope: ECon adds an edge and EJoin joins two ropes, so the
43	# two halves of a split join in one step.
44	type Edges is Data:
45	  ENil{}
46	  ECon{src: U32, dst: U32, w: F32, next: Edges}
47	  EJoin{a: Edges, b: Edges}
48	
49	# The options of ts2vg's graph classes: the visibility, how its edges point,
50	# their weight, the weights to keep (above `min_weight` and below
51	# `max_weight`; neg_inf() and inf() keep every weight), and how many points
52	# an edge may cross.
53	type Options is Data:
54	  Options{visibility: Visibility, directed: Direction, weighted: Weight, min_weight: F32, max_weight: F32, penetrable_limit: U32}
55	
56	# The `n` values `ys` of a time series, indexed from 0, and their x: the
57	# index, or `At{xs}` with the x of each value, increasing.
58	type Series is Type:
59	  Series{n: U32, ys: Array<F32>, xs: Xs}
60	
61	# A built graph: its `n` nodes, how its edges point, and its edges.
62	type Graph is Data:
63	  Graph{n: U32, directed: Direction, edges: Edges}
64	
65	# The graph as compressed rows: the neighbors of node i are
66	# `targets[offsets[i]]` up to `targets[offsets[i + 1]]`, increasing, and
67	# the weights of their edges sit at the same positions. An undirected
68	# graph's edge appears in the rows of both its ends, a directed one's in
69	# its source's row.
70	type Adjacency is Type:
71	  Adjacency{offsets: Array<U32>, targets: Array<U32>, weights: Array<F32>}
72	
73	# F32 constants, from divisions by zero: the JS lane builds an F32 from its
74	# bits with two typed arrays, too slow for a sweep that passes NaN at every
75	# point it visits.
76	
77	def nan() -> F32:
78	  F32.div(0.0, 0.0)
79	
80	def inf() -> F32:
81	  F32.div(1.0, 0.0)
82	
83	def neg_inf() -> F32:
84	  F32.div(F32.neg(1.0), 0.0)
85	
86	def is_nan(+a: F32) -> Bool:
87	  F32.is_ne(a, a)
88	
89	def pick(c: Bool, a: F32, b: F32) -> F32:
90	  match c:
91	    case True{}:
92	      a
93	    case False{}:
94	      b
95	
96	def pick_u32(c: Bool, a: U32, b: U32) -> U32:
97	  match c:
98	    case True{}:
99	      a
100	    case False{}:
101	      b
102	
103	# Weights: ts2vg's _weight_* functions. `s` is the slope the sweep measured,
104	# or NaN when it measured none.
105	
106	def sq_distance(+x1: F32, +x2: F32, +y1: F32, +y2: F32) -> F32:
107	  F32.add(F32.mul(F32.sub(x2, x1), F32.sub(x2, x1)), F32.mul(F32.sub(y2, y1), F32.sub(y2, y1)))
108	
109	def known_slope(+s: F32, +x1: F32, +x2: F32, +y1: F32, +y2: F32) -> F32:
110	  pick(is_nan(s), F32.div(F32.sub(y2, y1), F32.sub(x2, x1)), s)
111	
112	def weight(+k: Weight, +x1: F32, +x2: F32, +y1: F32, +y2: F32, +s: F32) -> F32:
113	  match k:
114	    case Unweighted{}:
115	      0.0
116	    case Distance{}:
117	      F32.sqrt(sq_distance(x1, x2, y1, y2))
118	    case SqDistance{}:
119	      sq_distance(x1, x2, y1, y2)
120	    case VDistance{}:
121	      F32.sub(y2, y1)
122	    case AbsVDistance{}:
123	      F32.abs(F32.sub(y2, y1))
124	    case HDistance{}:
125	      F32.sub(x2, x1)
126	    case AbsHDistance{}:
127	      F32.abs(F32.sub(x2, x1))
128	    case Slope{}:
129	      known_slope(s, x1, x2, y1, y2)
130	    case AbsSlope{}:
131	      F32.abs(known_slope(s, x1, x2, y1, y2))
132	    case Angle{}:
133	      F32.atan(known_slope(s, x1, x2, y1, y2))
134	    case AbsAngle{}:
135	      F32.atan(F32.abs(known_slope(s, x1, x2, y1, y2)))
136	    case NumPenetrations{}:
137	      0.0
138	
139	# Edges: ts2vg's add_edge, which drops an edge whose weight is at most lo or
140	# at least hi.
141	
142	def put(drop: Bool, i1: U32, i2: U32, w: F32, es: Edges) -> Edges:
143	  match drop:
144	    case True{}:
145	      es
146	    case False{}:
147	      ECon{i1, i2, w, es}
148	
149	def add_edge(+k: Weight, +lo: F32, +hi: F32, i1: U32, i2: U32, +x1: F32, +x2: F32, +y1: F32, +y2: F32, +s: F32, es: Edges) -> Edges:
150	  +w = weight(k, x1, x2, y1, y2, s)
151	  put(F32.is_le(w, lo) || F32.is_ge(w, hi), i1, i2, w, es)
152	
153	# A visible point right of the pivot a: the edge (a, b).
154	def add_right(seen: Bool, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, xa: F32, xb: F32, ya: F32, yb: F32, s: F32, es: Edges) -> Edges:
155	  match seen:
156	    case True{}:
157	      add_edge(k, lo, hi, ia, ib, xa, xb, ya, yb, s, es)
158	    case False{}:
159	      es
160	
161	# A visible point left of the pivot a: the edge (a, b) for top_to_bottom,
162	# else (b, a).
163	def add_left_dir(dir: Direction, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, xa: F32, xb: F32, ya: F32, yb: F32, s: F32, es: Edges) -> Edges:
164	  match dir:
165	    case TopToBottom{}:
166	      add_edge(k, lo, hi, ia, ib, xa, xb, ya, yb, s, es)
167	    case _:
168	      add_edge(k, lo, hi, ib, ia, xb, xa, yb, ya, s, es)
169	
170	def add_left(seen: Bool, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, xa: F32, xb: F32, ya: F32, yb: F32, s: F32, es: Edges) -> Edges:
171	  match seen:
172	    case True{}:
173	      add_left_dir(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, s, es)
174	    case False{}:
175	      es
176	
177	# Sweeps from the pivot a over the points of one side, nearest first.
178	
179	def abs_tol() -> F32:
180	  {0.00000000000001 : F32}
181	
182	def rel_tol() -> F32:
183	  {0.00000000000001 : F32}
184	
185	# max(|xa|, |ya|), which a sweep computes once.
186	def scale(+xa: F32, +ya: F32) -> F32:
187	  F32.max(F32.abs(xa), F32.abs(ya))
188	
189	# ts2vg's tolerance for the points a and b.
190	# F32.max only picks one of its words, so grouping the four maxima this way
191	# gives ts2vg's word.
192	def tol(+scale_a: F32, +xb: F32, +yb: F32) -> F32:
193	  F32.max(abs_tol(), F32.mul(rel_tol(), F32.max(scale_a, F32.max(F32.abs(xb), F32.abs(yb)))))
194	
195	def greater(+a: F32, +b: F32, +t: F32) -> Bool:
196	  F32.is_gt(F32.sub(a, b), t)
197	
198	# greater(s, best, tol(scale_a, xb, yb)), with the tolerance computed only
199	# when s is above best: the tolerance is positive, so otherwise it fails
200	# anyway.
201	def clears(up: Bool, +s: F32, +best: F32, +scale_a: F32, +xb: F32, +yb: F32) -> Bool:
202	  match up:
203	    case True{}:
204	      greater(s, best, tol(scale_a, xb, yb))
205	    case False{}:
206	      False{}
207	
208	# Divide and conquer: ts2vg's _compute_graph takes an interval's first
209	# highest point as its pivot, sweeps the interval's points left and right of
210	# it, and recurses on both sides. Every point is the pivot of one interval,
211	# so once each point knows how many points its interval holds on each side,
212	# the sweeps no longer depend on each other.
213	
214	# The stack holds the points no later point has beaten yet. A point pops the
215	# entries it beats; the entry left on top bounds its interval on the left,
216	# and the point that pops an entry bounds that entry's interval on the right.
217	type Stack is Data:
218	  KNil{}
219	  KCon{i: U32, key: F32, next: Stack}
220	
221	# What a later point must beat to pop y. An equal point beats nothing, so an
222	# interval's first highest point is its pivot. A NaN beats nothing either,
223	# and ts2vg's argmax takes it only as an interval's first point, so a NaN at
224	# i is the pivot of the right side of i - 1: it takes the key of the entry
225	# below it and leaves the stack with that entry.
226	def key(+y: F32, +below: F32) -> F32:
227	  pick(is_nan(y), below, y)
228	
229	# The work only balances the split, so F32 rounding does no harm.
230	type Popped is Type:
231	  Popped{first: U32, below: F32, st: Stack, rights: Array<U32>, visits: F32}
232	
233	def pops(rest: Stack, beat: Bool, +y: F32, +i: U32, +j: U32, +kj: F32, rights: Array<U32>, +visits: F32) -> Popped:
234	  match rest beat:
235	    case KCon{+j2, +k2, rest2} True{}:
236	      +r = (i - 1 - j : U32)
237	      pops(rest2, F32.is_gt(y, k2), y, i, j2, k2, Array.set(U32, rights, j, r), F32.add(visits, U32.to_f32(r)))
238	    case KNil{} True{}:
239	      +r = (i - 1 - j : U32)
240	      Popped{0, inf(), KNil{}, Array.set(U32, rights, j, r), F32.add(visits, U32.to_f32(r))}
241	    case _ False{}:
242	      Popped{(j + 1 : U32), kj, KCon{j, kj, rest}, rights, visits}
243	
244	def flush(st: Stack, +n: U32, rights: Array<U32>, +visits: F32) -> Array<U32> & F32:
245	  match st:
246	    case KCon{+j, k, rest}:
247	      +r = (n - 1 - j : U32)
248	      flush(rest, n, Array.set(U32, rights, j, r), F32.add(visits, U32.to_f32(r)))
249	    case KNil{}:
250	      (rights, visits)
251	
252	# lefts and rights count the points each point's left and right sweeps visit.
253	type Cols is Type:
254	  Cols{ys: Array<F32>, lefts: Array<U32>, rights: Array<U32>}
255	
256	type Arrived is Type:
257	  OpensWithNaN{rights: Array<U32>, visits: F32}
258	  PopsAll{rights: Array<U32>, visits: F32}
259	  Bounded{first: U32, below: F32, st: Stack, rights: Array<U32>, visits: F32}
260	
261	def lpops(rest: Stack, beat: Bool, +y: F32, +i: U32, +j: U32, +kj: F32, rights: Array<U32>, +visits: F32) -> Arrived:
262	  match rest beat:
263	    case KCon{+j2, +k2, rest2} True{}:
264	      +r = (i - 1 - j : U32)
265	      lpops(rest2, F32.is_gt(y, k2), y, i, j2, k2, Array.set(U32, rights, j, r), F32.add(visits, U32.to_f32(r)))
266	    case KNil{} True{}:
267	      +r = (i - 1 - j : U32)
268	      PopsAll{Array.set(U32, rights, j, r), F32.add(visits, U32.to_f32(r))}
269	    case _ False{}:
270	      Bounded{(j + 1 : U32), kj, KCon{j, kj, rest}, rights, visits}
271	
272	def arrive(st: Stack, nan: Bool, +y: F32, +i: U32, rights: Array<U32>, +visits: F32) -> Arrived:
273	  match st nan:
274	    case KNil{} True{}:
275	      OpensWithNaN{rights, visits}
276	    case KNil{} False{}:
277	      PopsAll{rights, visits}
278	    case KCon{+j, +kj, rest} True{}:
279	      Bounded{(j + 1 : U32), kj, KCon{j, kj, rest}, rights, visits}
280	    case KCon{+j, +kj, rest} False{}:
281	      lpops(rest, F32.is_gt(y, kj), y, i, j, kj, rights, visits)
282	
283	# Latest first.
284	type Unbounded is Data:
285	  UNil{}
286	  UCon{i: U32, y: F32, next: Unbounded}
287	
288	# A block of points runs the stack on its own points. Its points that pop
289	# its whole stack, and the NaNs that open it, need the points before it.
290	type Block is Data:
291	  Block{start: U32, opening_nans: U32, unbounded: Unbounded, st: Stack, visits: F32}
292	
293	type Local is Type:
294	  Local{b: Block, c: Cols}
295	
296	def local(count: Nat, +i: U32, +y: F32, a: Arrived, yr: Array<F32> & F32, lefts: Array<U32>, +start: U32, +opening_nans: U32, unbounded: Unbounded) -> Local:
297	  match count a:
298	    case 1n+c OpensWithNaN{rights, +visits}:
299	      (ys, +y2) = yr
300	      +i2 = (i + 1 : U32)
301	      local(c, i2, y2, arrive(KNil{}, is_nan(y2), y2, i2, rights, F32.add(visits, 1.0)), Array.get(F32, ys, (i2 + 1 : U32)), Array.set(U32, lefts, i, 0), start, (opening_nans + 1 : U32), unbounded)
302	    case 1n+c PopsAll{rights, +visits}:
303	      (ys, +y2) = yr
304	      +i2 = (i + 1 : U32)
305	      local(c, i2, y2, arrive(KCon{i, y, KNil{}}, is_nan(y2), y2, i2, rights, visits), Array.get(F32, ys, (i2 + 1 : U32)), lefts, start, opening_nans, UCon{i, y, unbounded})
306	    case 1n+c Bounded{+first, +below, st, rights, +visits}:
307	      (ys, +y2) = yr
308	      +i2 = (i + 1 : U32)
309	      +l = (i - first : U32)
310	      local(c, i2, y2, arrive(KCon{i, key(y, below), st}, is_nan(y2), y2, i2, rights, F32.add(visits, U32.to_f32((l + 1 : U32)))), Array.get(F32, ys, (i2 + 1 : U32)), Array.set(U32, lefts, i, l), start, opening_nans, unbounded)
311	    case 0n OpensWithNaN{rights, +visits}:
312	      (ys, y2) = yr
313	      Local{Block{start, (opening_nans + 1 : U32), unbounded, KNil{}, F32.add(visits, 1.0)}, Cols{ys, Array.set(U32, lefts, i, 0), rights}}
314	    case 0n PopsAll{rights, +visits}:
315	      (ys, y2) = yr
316	      Local{Block{start, opening_nans, UCon{i, y, unbounded}, KCon{i, y, KNil{}}, visits}, Cols{ys, lefts, rights}}
317	    case 0n Bounded{+first, +below, st, rights, +visits}:
318	      (ys, y2) = yr
319	      +l = (i - first : U32)
320	      Local{Block{start, opening_nans, unbounded, KCon{i, key(y, below), st}, F32.add(visits, U32.to_f32((l + 1 : U32)))}, Cols{ys, Array.set(U32, lefts, i, l), rights}}
321	
322	def block_pass(count: Nat, +s: U32, yr: Array<F32> & F32, lefts: Array<U32>, rights: Array<U32>) -> Local:
323	  match count:
324	    case 0n:
325	      (ys, y) = yr
326	      Local{Block{s, 0, UNil{}, KNil{}, 0.0}, Cols{ys, lefts, rights}}
327	    case 1n+c:
328	      (ys, +y) = yr
329	      local(c, s, y, arrive(KNil{}, is_nan(y), y, s, rights, 0.0), Array.get(F32, ys, (s + 1 : U32)), lefts, s, 0, UNil{})
330	
331	type Merged is Type:
332	  Merged{st: Stack, lefts: Array<U32>, rights: Array<U32>, visits: F32}
333	
334	def push_with_key(count: Nat, +i: U32, +k: F32, st: Stack) -> Stack:
335	  match count:
336	    case 0n:
337	      st
338	    case 1n+c:
339	      push_with_key(c, (i + 1 : U32), k, KCon{i, k, st})
340	
341	def push_opening_nans(count: Nat, +i: U32, st: Stack) -> Stack:
342	  match st:
343	    case KCon{+j, +kj, rest}:
344	      push_with_key(count, i, kj, KCon{j, kj, rest})
345	    case KNil{}:
346	      push_with_key(count, i, inf(), KNil{})
347	
348	def pops_from(st: Stack, +y: F32, +i: U32, rights: Array<U32>, +visits: F32) -> Popped:
349	  match st:
350	    case KCon{+j, +kj, rest}:
351	      pops(rest, F32.is_gt(y, kj), y, i, j, kj, rights, visits)
352	    case KNil{}:
353	      Popped{0, inf(), KNil{}, rights, visits}
354	
355	def set_left(+i: U32, lefts: Array<U32>, p: Popped) -> Merged:
356	  Popped{+first, below, st, rights, +visits} = p
357	  +l = (i - first : U32)
358	  Merged{st, Array.set(U32, lefts, i, l), rights, F32.add(visits, U32.to_f32((l + 1 : U32)))}
359	
360	def in_order(u: Unbounded, acc: Unbounded) -> Unbounded:
361	  match u:
362	    case UNil{}:
363	      acc
364	    case UCon{i, y, next}:
365	      in_order(next, UCon{i, y, acc})
366	
367	def bound_unbounded(u: Unbounded, m: Merged) -> Merged:
368	  match u:
369	    case UNil{}:
370	      m
371	    case UCon{+i, +y, next}:
372	      Merged{st, lefts, rights, +visits} = m
373	      bound_unbounded(next, set_left(i, lefts, pops_from(st, y, i, rights, visits)))
374	
375	def onto(st: Stack, acc: Stack) -> Stack:
376	  match st:
377	    case KNil{}:
378	      acc
379	    case KCon{j, k, rest}:
380	      onto(rest, KCon{j, k, acc})
381	
382	def push_block_stack(st: Stack, m: Merged) -> Merged:
383	  Merged{g, lefts, rights, visits} = m
384	  Merged{onto(onto(st, KNil{}), g), lefts, rights, visits}
385	
386	# The merge replays the stack over every point before the block. The points
387	# the block's own stack could not bound pop it in order; the others sit on
388	# the block's stack, which goes on top.
389	def merge_block(b: Block, m: Merged) -> Merged:
390	  Block{+start, +opening_nans, unbounded, st, +visits} = b
391	  Merged{g, lefts, rights, +v} = m
392	  push_block_stack(st, bound_unbounded(in_order(unbounded, UNil{}), Merged{push_opening_nans(U32.to_nat(opening_nans), start, g), lefts, rights, F32.add(v, visits)}))
393	
394	# Bend matches only parameters, so a read's result must arrive as one: a
395	# sweep reads y and x at `ahead`, where the next sweep starts, and hands them
396	# over. The sweeps read x through `x_at`, so x as the index costs no read.
397	type Walked is Type:
398	  Walked{yr: Array<F32> & F32, xr: Array<F32> & F32, es: Edges}
399	
400	def at_index(xs: Array<F32>, +j: U32) -> Array<F32> & F32:
401	  (xs, U32.to_f32(j))
402	
403	def natural_right(~x_at: Array<F32> -> U32 -> Array<F32> & F32, count: Nat, +j: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +scale_a: F32, +best: F32, es: Edges, +ahead: U32) -> Walked:
404	  match count:
405	    case 0n:
406	      (ys, yb) = yr
407	      (xs, xb) = xr
408	      Walked{Array.get(F32, ys, ahead), Array.get(F32, xs, ahead), es}
409	    case 1n+c:
410	      (ys, +yb) = yr
411	      (xs, +xb) = xr
412	      +s = F32.div(F32.sub(yb, ya), F32.sub(xb, xa))
413	      +seen = clears(F32.is_gt(s, best), s, best, scale_a, xb, yb)
414	      +j2 = (j + 1 : U32)
415	      natural_right(~x_at, c, j2, Array.get(F32, ys, j2), x_at(xs, j2), k, lo, hi, ia, xa, ya, scale_a, pick(seen, s, best), add_right(seen, k, lo, hi, ia, j, xa, xb, ya, yb, s, es), ahead)
416	
417	def natural_left(~x_at: Array<F32> -> U32 -> Array<F32> & F32, count: Nat, +j: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +scale_a: F32, +best: F32, es: Edges, +ahead: U32) -> Walked:
418	  match count:
419	    case 0n:
420	      (ys, yb) = yr
421	      (xs, xb) = xr
422	      Walked{Array.get(F32, ys, ahead), Array.get(F32, xs, ahead), es}
423	    case 1n+c:
424	      (ys, +yb) = yr
425	      (xs, +xb) = xr
426	      +s = F32.div(F32.sub(yb, ya), F32.neg(F32.sub(xb, xa)))
427	      +seen = clears(F32.is_gt(s, best), s, best, scale_a, xb, yb)
428	      +j2 = (j - 1 : U32)
429	      natural_left(~x_at, c, j2, Array.get(F32, ys, j2), x_at(xs, j2), dir, k, lo, hi, ia, xa, ya, scale_a, pick(seen, s, best), add_left(seen, dir, k, lo, hi, ia, j, xa, xb, ya, yb, F32.neg(s), es), ahead)
430	
431	def horizontal_right(~x_at: Array<F32> -> U32 -> Array<F32> & F32, count: Nat, +j: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +best: F32, es: Edges, +ahead: U32) -> Walked:
432	  match count:
433	    case 0n:
434	      (ys, yb) = yr
435	      (xs, xb) = xr
436	      Walked{Array.get(F32, ys, ahead), Array.get(F32, xs, ahead), es}
437	    case 1n+c:
438	      (ys, +yb) = yr
439	      (xs, xb) = xr
440	      +seen = F32.is_gt(yb, best)
441	      +j2 = (j + 1 : U32)
442	      horizontal_right(~x_at, c, j2, Array.get(F32, ys, j2), x_at(xs, j2), k, lo, hi, ia, xa, ya, pick(seen, yb, best), add_right(seen, k, lo, hi, ia, j, xa, xb, ya, yb, nan(), es), ahead)
443	
444	def horizontal_left(~x_at: Array<F32> -> U32 -> Array<F32> & F32, count: Nat, +j: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +best: F32, es: Edges, +ahead: U32) -> Walked:
445	  match count:
446	    case 0n:
447	      (ys, yb) = yr
448	      (xs, xb) = xr
449	      Walked{Array.get(F32, ys, ahead), Array.get(F32, xs, ahead), es}
450	    case 1n+c:
451	      (ys, +yb) = yr
452	      (xs, xb) = xr
453	      +seen = F32.is_gt(yb, best)
454	      +j2 = (j - 1 : U32)
455	      horizontal_left(~x_at, c, j2, Array.get(F32, ys, j2), x_at(xs, j2), dir, k, lo, hi, ia, xa, ya, pick(seen, yb, best), add_left(seen, dir, k, lo, hi, ia, j, xa, xb, ya, yb, nan(), es), ahead)
456	
457	# Each sweep, with x as the index (`index`) or read from `xs`.
458	def right_sweep(g: Visibility, index: Bool, count: Nat, yr: Array<F32> & F32, xs: Array<F32>, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, es: Edges) -> Walked:
459	  match g index:
460	    case Natural{} True{}:
461	      natural_right(~(a => +j => at_index(a, j)), count, (ia + 1 : U32), yr, at_index(xs, (ia + 1 : U32)), k, lo, hi, ia, xa, ya, scale(xa, ya), neg_inf(), es, (ia - 1 : U32))
462	    case Natural{} False{}:
463	      natural_right(~(a => +j => Array.get(F32, a, j)), count, (ia + 1 : U32), yr, Array.get(F32, xs, (ia + 1 : U32)), k, lo, hi, ia, xa, ya, scale(xa, ya), neg_inf(), es, (ia - 1 : U32))
464	    case Horizontal{} True{}:
465	      horizontal_right(~(a => +j => at_index(a, j)), count, (ia + 1 : U32), yr, at_index(xs, (ia + 1 : U32)), k, lo, hi, ia, xa, ya, neg_inf(), es, (ia - 1 : U32))
466	    case Horizontal{} False{}:
467	      horizontal_right(~(a => +j => Array.get(F32, a, j)), count, (ia + 1 : U32), yr, Array.get(F32, xs, (ia + 1 : U32)), k, lo, hi, ia, xa, ya, neg_inf(), es, (ia - 1 : U32))
468	
469	def left_sweep(g: Visibility, index: Bool, count: Nat, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, w: Walked) -> Walked:
470	  match g index:
471	    case Natural{} True{}:
472	      Walked{yr, xr, es} = w
473	      (xs, xv) = xr
474	      natural_left(~(a => +j => at_index(a, j)), count, (ia - 1 : U32), yr, at_index(xs, (ia - 1 : U32)), dir, k, lo, hi, ia, xa, ya, scale(xa, ya), neg_inf(), es, (ia + 1 : U32))
475	    case Natural{} False{}:
476	      Walked{yr, xr, es} = w
477	      (xs, xv) = xr
478	      natural_left(~(a => +j => Array.get(F32, a, j)), count, (ia - 1 : U32), yr, Array.get(F32, xs, (ia - 1 : U32)), dir, k, lo, hi, ia, xa, ya, scale(xa, ya), neg_inf(), es, (ia + 1 : U32))
479	    case Horizontal{} True{}:
480	      Walked{yr, xr, es} = w
481	      (xs, xv) = xr
482	      horizontal_left(~(a => +j => at_index(a, j)), count, (ia - 1 : U32), yr, at_index(xs, (ia - 1 : U32)), dir, k, lo, hi, ia, xa, ya, neg_inf(), es, (ia + 1 : U32))
483	    case Horizontal{} False{}:
484	      Walked{yr, xr, es} = w
485	      (xs, xv) = xr
486	      horizontal_left(~(a => +j => Array.get(F32, a, j)), count, (ia - 1 : U32), yr, Array.get(F32, xs, (ia - 1 : U32)), dir, k, lo, hi, ia, xa, ya, neg_inf(), es, (ia + 1 : U32))
487	
488	type Out is Type:
489	  Out{cols: Cols, xs: Array<F32>, es: Edges}
490	
491	def points(count: Nat, +i: U32, lr: Array<U32> & U32, rr: Array<U32> & U32, w: Walked, +index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32) -> Out:
492	  match count:
493	    case 0n:
494	      (lefts, l) = lr
495	      (rights, r) = rr
496	      Walked{yr, xr, es} = w
497	      (ys, y) = yr
498	      (xs, x) = xr
499	      Out{Cols{ys, lefts, rights}, xs, es}
500	    case 1n+c:
501	      (lefts, +l) = lr
502	      (rights, +r) = rr
503	      Walked{yr, xr, es} = w
504	      (ys, +y) = yr
505	      (xs, +xv) = xr
506	      +x = pick(index, U32.to_f32(i), xv)
507	      +i2 = (i + 1 : U32)
508	      points(c, i2, Array.get(U32, lefts, i2), Array.get(U32, rights, i2), left_sweep(g, index, U32.to_nat(l), dir, k, lo, hi, i, x, y, right_sweep(g, index, U32.to_nat(r), Array.get(F32, ys, i2), xs, k, lo, hi, i, x, y, es)), index, g, dir, k, lo, hi)
509	
510	def edges_of(o: Out) -> Edges:
511	  Out{cols, xs, es} = o
512	  es
513	
514	def swept(+index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +n: U32, ys: Array<F32>, xs: Array<F32>, lefts: Array<U32>, r: Array<U32> & F32) -> Edges:
515	  (rights, visits) = r
516	  edges_of(points(U32.to_nat(n), 0, Array.get(U32, lefts, 0), Array.get(U32, rights, 0), Walked{Array.get(F32, ys, 0), Array.get(F32, xs, 0), ENil{}}, index, g, dir, k, lo, hi))
517	
518	def bounded(+index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +n: U32, ys: Array<F32>, xs: Array<F32>, m: Merged) -> Edges:
519	  Merged{st, lefts, rights, +visits} = m
520	  swept(index, g, dir, k, lo, hi, n, ys, xs, lefts, flush(st, n, rights, visits))
521	
522	def stacked(+index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +n: U32, xs: Array<F32>, l: Local) -> Edges:
523	  Local{b, c} = l
524	  Cols{ys, lefts, rights} = c
525	  bounded(index, g, dir, k, lo, hi, n, ys, xs, merge_block(b, Merged{KNil{}, lefts, rights, 0.0}))
526	
527	# An array with room for `n` slots, indexed from 0.
528	def slots(+n: U32) -> Array<U32>:
529	  Array.new(U32, Nat.add(1n, U32.log2(n)), 0)
530	
531	# A stand-in for `xs` when x is the index, which no sweep reads.
532	def no_xs() -> Array<F32>:
533	  Array.new(F32, 0n, 0.0)
534	
535	def built(+index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, ys: Array<F32>, xs: Array<F32>, +n: U32) -> Edges:
536	  stacked(index, g, dir, k, lo, hi, n, xs, block_pass(U32.to_nat(n), 0, Array.get(F32, ys, 0), slots(n), slots(n)))
537	
538	# The graph of the `n` values `ys` at `xs`, on one thread, as one block.
539	# `lo` and `hi` bound the weights, as ts2vg's min_weight and max_weight;
540	# pass neg_inf() and inf() for no bound. The x of `Index{}` is the index as
541	# an F32, exact for up to 2^24 values; the x of `At{xs}` must increase.
542	# parallel.bend builds the same graph on every thread, and relies on
543	# Base's @unsafe Array.fork.
544	def visibility(+g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, ys: Array<F32>, xs: Xs, +n: U32) -> Edges:
545	  match xs:
546	    case Index{}:
547	      built(True{}, g, dir, k, lo, hi, ys, no_xs(), n)
548	    case At{a}:
549	      built(False{}, g, dir, k, lo, hi, ys, a, n)
550	
551	# Degrees: each edge adds one to the degree of both of its ends, as ts2vg's
552	# degrees_in + degrees_out.
553	
554	def added(r: Array<U32> & U32) -> Array<U32>:
555	  (d, old) = r
556	  d
557	
558	def atomic_bump(d: Array<U32>, +i: U32) -> Array<U32>:
559	  added(Array.atomic.add(d, i, 1))
560	
561	# The degrees of the nodes of `es`, added into `d`; parallel.bend counts
562	# them on every thread.
563	def degrees.into(es: Edges, d: Array<U32>) -> Array<U32>:
564	  match es:
565	    case ENil{}:
566	      d
567	    case ECon{s, t, w, next}:
568	      degrees.into(next, atomic_bump(atomic_bump(d, s), t))
569	    case EJoin{a, b}:
570	      degrees.into(b, degrees.into(a, d))
571	
572	# Limited penetrable visibility graphs, as ts2vg 1.2.4's _natural_penetrable
573	# and _horizontal_penetrable build them, in F32: every point sweeps the
574	# points right of it, and an edge may cross up to `limit` points. The sweeps
575	# are independent, so the points split into parallel ranges of equal work.
576	
577	# The points, in a list the parallel ranges share.
578	type Pts is Data:
579	  PNil{}
580	  PCon{i: U32, x: F32, y: F32, next: Pts}
581	
582	# The points of `ys` and `xs` from index `i` down to 0, onto `acc`; `yr`
583	# and `xr` hold y and x at `i`, and x is the index if `index` says so. The
584	# last read, at index -1, wraps around.
585	def penetrable.points(count: Nat, +i: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +index: Bool, acc: Pts) -> Pts:
586	  match count:
587	    case 0n:
588	      (ys, y) = yr
589	      (xs, x) = xr
590	      acc
591	    case 1n+c:
592	      (ys, +y) = yr
593	      (xs, +x) = xr
594	      +i2 = (i - 1 : U32)
595	      penetrable.points(c, i2, Array.get(F32, ys, i2), Array.get(F32, xs, i2), index, PCon{i, pick(index, U32.to_f32(i), x), y, acc})
596	
597	# The `limit + 1` largest slopes (or heights) a sweep has accepted, as ts2vg's
598	# max_slopes; the sweep accepts a point above the smallest of them.
599	type Tops is Data:
600	  TNil{}
601	  TCon{v: F32, next: Tops}
602	
603	# The position and value of the first smallest top.
604	type Low is Data:
605	  Low{k: Nat, v: F32}
606	
607	def tops_new(n: Nat) -> Tops:
608	  match n:
609	    case 0n:
610	      TCon{neg_inf(), TNil{}}
611	    case 1n+p:
612	      TCon{neg_inf(), tops_new(p)}
613	
614	def pick_low(c: Bool, a: Low, b: Low) -> Low:
615	  match c:
616	    case True{}:
617	      a
618	    case False{}:
619	      b
620	
621	def low_from(ts: Tops, +k: Nat, +best: Low) -> Low:
622	  match ts:
623	    case TNil{}:
624	      best
625	    case TCon{+v, next}:
626	      Low{bk, +bv} = best
627	      low_from(next, 1n+k, pick_low(F32.is_lt(v, bv), Low{k, v}, Low{bk, bv}))
628	
629	def lowest(ts: Tops) -> Low:
630	  low_from(ts, 0n, Low{0n, inf()})
631	
632	def Low.v(+l: Low) -> F32:
633	  Low{k, v} = l
634	  v
635	
636	def set_at(k: Nat, ts: Tops, v: F32) -> Tops:
637	  match k ts:
638	    case 0n TCon{old, next}:
639	      TCon{v, next}
640	    case 1n+p TCon{old, next}:
641	      TCon{old, set_at(p, next, v)}
642	    case _ TNil{}:
643	      TNil{}
644	
645	# How many tops a natural edge with slope `s` does not clear, as ts2vg's
646	# num_penetrations weight.
647	def blocked_natural(ts: Tops, +s: F32, +t: F32, +n: F32) -> F32:
648	  match ts:
649	    case TNil{}:
650	      n
651	    case TCon{v, next}:
652	      blocked_natural(next, s, t, pick(greater(s, v, t), n, F32.add(n, 1.0)))
653	
654	def blocked_horizontal(ts: Tops, +ya: F32, +yb: F32, +n: F32) -> F32:
655	  match ts:
656	    case TNil{}:
657	      n
658	    case TCon{+v, next}:
659	      blocked_horizontal(next, ya, yb, pick(F32.is_le(ya, v) || F32.is_le(yb, v), F32.add(n, 1.0), n))
660	
661	# ts2vg's add_edge with a known weight: NaN means the weight function
662	# computes it.
663	def add_known(+k: Weight, +lo: F32, +hi: F32, i1: U32, i2: U32, +x1: F32, +x2: F32, +y1: F32, +y2: F32, +s: F32, +known: F32, es: Edges) -> Edges:
664	  +w = pick(is_nan(known), weight(k, x1, x2, y1, y2, s), known)
665	  put(F32.is_le(w, lo) || F32.is_ge(w, hi), i1, i2, w, es)
666	
667	# The edge from a to b, or from b to a for top_to_bottom when b is higher.
668	def add_oriented(down: Bool, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, xa: F32, xb: F32, ya: F32, yb: F32, s: F32, known: F32, es: Edges) -> Edges:
669	  match down:
670	    case True{}:
671	      add_known(k, lo, hi, ib, ia, xb, xa, yb, ya, s, known, es)
672	    case False{}:
673	      add_known(k, lo, hi, ia, ib, xa, xb, ya, yb, s, known, es)
674	
675	def is_down(dir: Direction, +ya: F32, +yb: F32) -> Bool:
676	  match dir:
677	    case TopToBottom{}:
678	      F32.is_gt(yb, ya)
679	    case _:
680	      False{}
681	
682	def is_count(k: Weight) -> Bool:
683	  match k:
684	    case NumPenetrations{}:
685	      True{}
686	    case _:
687	      False{}
688	
689	def natural_edge(+dir: Direction, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, +xa: F32, +xb: F32, +ya: F32, +yb: F32, +s: F32, +scale_a: F32, ts: Tops, es: Edges) -> Edges:
690	  +known = pick(is_count(k), blocked_natural(ts, s, tol(scale_a, xb, yb), 0.0), nan())
691	  add_oriented(is_down(dir, ya, yb), k, lo, hi, ia, ib, xa, xb, ya, yb, s, known, es)
692	
693	def horizontal_edge(+dir: Direction, +k: Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, +xa: F32, xb: F32, +ya: F32, +yb: F32, ts: Tops, es: Edges) -> Edges:
694	  +known = pick(is_count(k), blocked_horizontal(ts, ya, yb, 0.0), nan())
695	  add_oriented(is_down(dir, ya, yb), k, lo, hi, ia, ib, xa, xb, ya, yb, nan(), known, es)
696	
697	# The sweep of point a over the points right of it: ts2vg's inner loop. A
698	# point b the sweep sees replaces the smallest top and adds its edge at the
699	# next step, when `up` says so, so a step after an unseen point calls no
700	# def that recurses.
701	def natural_sweep(pts: Pts, up: Bool, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +scale_a: F32, +ts: Tops, low: Low, +ib: U32, +xb: F32, +yb: F32, +sb: F32, es: Edges) -> Edges:
702	  match pts up:
703	    case PNil{} False{}:
704	      es
705	    case PNil{} True{}:
706	      natural_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, sb, scale_a, ts, es)
707	    case PCon{+ic, +xc, +yc, next} False{}:
708	      Low{+lk, +lv} = low
709	      +s = F32.div(F32.sub(yc, ya), F32.sub(xc, xa))
710	      natural_sweep(next, clears(F32.is_gt(s, lv), s, lv, scale_a, xc, yc), dir, k, lo, hi, ia, xa, ya, scale_a, ts, Low{lk, lv}, ic, xc, yc, s, es)
711	    case PCon{+ic, +xc, +yc, next} True{}:
712	      Low{+lk, lv} = low
713	      +after = set_at(lk, ts, sb)
714	      +now = lowest(after)
715	      +s = F32.div(F32.sub(yc, ya), F32.sub(xc, xa))
716	      natural_sweep(next, clears(F32.is_gt(s, Low.v(now)), s, Low.v(now), scale_a, xc, yc), dir, k, lo, hi, ia, xa, ya, scale_a, after, now, ic, xc, yc, s, natural_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, sb, scale_a, ts, es))
717	
718	# The same for the horizontal sweep, which stops once the smallest top is
719	# above a: no later point can be seen.
720	def horizontal_sweep(pts: Pts, up: Bool, stop: Bool, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +ts: Tops, low: Low, +ib: U32, +xb: F32, +yb: F32, es: Edges) -> Edges:
721	  match pts up stop:
722	    case _ False{} True{}:
723	      es
724	    case PNil{} False{} False{}:
725	      es
726	    case PNil{} True{} _:
727	      horizontal_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, ts, es)
728	    case PCon{+ic, xc, +yc, next} False{} False{}:
729	      Low{+lk, +lv} = low
730	      horizontal_sweep(next, F32.is_gt(ya, lv) && F32.is_gt(yc, lv), False{}, dir, k, lo, hi, ia, xa, ya, ts, Low{lk, lv}, ic, xc, yc, es)
731	    case PCon{+ic, xc, +yc, next} True{} _:
732	      Low{+lk, lv} = low
733	      +after = set_at(lk, ts, yb)
734	      +now = lowest(after)
735	      horizontal_sweep(next, F32.is_gt(ya, Low.v(now)) && F32.is_gt(yc, Low.v(now)), F32.is_gt(Low.v(now), ya), dir, k, lo, hi, ia, xa, ya, after, now, ic, xc, yc, horizontal_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, ts, es))
736	
737	def sweep(g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: Nat, ia: U32, +xa: F32, +ya: F32, pts: Pts, es: Edges) -> Edges:
738	  match g:
739	    case Natural{}:
740	      natural_sweep(pts, False{}, dir, k, lo, hi, ia, xa, ya, scale(xa, ya), tops_new(limit), Low{0n, neg_inf()}, 0, 0.0, 0.0, 0.0, es)
741	    case Horizontal{}:
742	      horizontal_sweep(pts, False{}, False{}, dir, k, lo, hi, ia, xa, ya, tops_new(limit), Low{0n, neg_inf()}, 0, 0.0, 0.0, es)
743	
744	# The sweeps of the first `count` points of `pts`, one rope each: ts2vg's
745	# outer loop.
746	def sweeps(pts: Pts, count: Nat, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: Nat) -> Edges:
747	  match pts count:
748	    case PCon{i, x, y, +next} 1n+c:
749	      EJoin{sweep(g, dir, k, lo, hi, limit, i, x, y, next, ENil{}), sweeps(next, c, g, dir, k, lo, hi, limit)}
750	    case _ _:
751	      ENil{}
752	
753	def drop(n: Nat, pts: Pts) -> Pts:
754	  match n pts:
755	    case 1n+p PCon{i, x, y, next}:
756	      drop(p, next)
757	    case _ _:
758	      pts
759	
760	# How many of the first `count` of `len` points do half the sweeps' work. A
761	# natural sweep walks every point right of its start, so the first c starts
762	# walk c (2 len - 1 - c) / 2 points; a horizontal sweep stops early, so
763	# halving the starts is close enough.
764	def first_half(g: Visibility, +len: U32, +count: U32) -> U32:
765	  match g:
766	    case Natural{}:
767	      +b = F32.sub(F32.mul(2.0, U32.to_f32(len)), 1.0)
768	      +c = U32.to_f32(count)
769	      U32.min(count, F32.to_u32(F32.div(F32.sub(b, F32.sqrt(F32.sub(F32.mul(b, b), F32.mul(2.0, F32.mul(c, F32.sub(b, c)))))), 2.0)))
770	    case Horizontal{}:
771	      U32.div(count, 2)
772	
773	# The sweeps of the first `count` points of `pts`, as two halves of equal
774	# work in parallel, `depth` times over. `len` and `cu` are `pts`'s length
775	# and `count` as U32; they only place the cut.
776	def split(depth: Nat, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: Nat, +pts: Pts, +len: U32, +cu: U32, +count: Nat) -> Edges:
777	  match depth:
778	    case 0n:
779	      sweeps(pts, count, g, dir, k, lo, hi, limit)
780	    case 1n++d:
781	      +c1u = first_half(g, len, cu)
782	      +c1 = Nat.min(count, U32.to_nat(c1u))
783	      a b = split(d, g, dir, k, lo, hi, limit, pts, len, c1u, c1) split(d, g, dir, k, lo, hi, limit, drop(c1, pts), (len - c1u : U32), (cu - c1u : U32), Nat.sub(count, c1))
784	      EJoin{a, b}
785	
786	def penetrable.built(+index: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: U32, ys: Array<F32>, xs: Array<F32>, +n: U32) -> Edges:
787	  split(7n, g, dir, k, lo, hi, U32.to_nat(limit), penetrable.points(U32.to_nat(n), (n - 1 : U32), Array.get(F32, ys, (n - 1 : U32)), Array.get(F32, xs, (n - 1 : U32)), index, PNil{}), n, n, U32.to_nat(n))
788	
789	# The limited penetrable graph of the `n` values `ys` at `xs`, crossing at
790	# most `limit` points. The x of `Index{}` is the index as an F32, exact
791	# for up to 2^24 values; the x of `At{xs}` must increase.
792	def penetrable(+g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: U32, ys: Array<F32>, xs: Xs, +n: U32) -> Edges:
793	  match xs:
794	    case Index{}:
795	      penetrable.built(True{}, g, dir, k, lo, hi, limit, ys, no_xs(), n)
796	    case At{a}:
797	      penetrable.built(False{}, g, dir, k, lo, hi, limit, ys, a, n)
798	
799	# The library
800	
801	# ts2vg's defaults for the visibility `v`: undirected, unweighted, every
802	# weight kept, no point crossed.
803	def options(v: Visibility) -> Options:
804	  Options{v, Undirected{}, Unweighted{}, neg_inf(), inf(), 0}
805	
806	def length(vs: List<&2, F32>, +n: U32) -> U32:
807	  match vs:
808	    case Nil{}:
809	      n
810	    case Con{v, rest}:
811	      length(rest, (n + 1 : U32))
812	
813	def fill(vs: List<&2, F32>, +i: U32, a: Array<F32>) -> Array<F32>:
814	  match vs:
815	    case Nil{}:
816	      a
817	    case Con{v, rest}:
818	      fill(rest, (i + 1 : U32), Array.set(F32, a, i, v))
819	
820	def column(+n: U32, vs: List<&2, F32>) -> Array<F32>:
821	  fill(vs, 0, Array.new(F32, Nat.add(1n, U32.log2(n)), 0.0))
822	
823	def counted(+n: U32, ys: List<&2, F32>) -> Series:
824	  Series{n, column(n, ys), Index{}}
825	
826	# The series of the values `ys`, each at its index.
827	def series(+ys: List<&2, F32>) -> Series:
828	  counted(length(ys, 0), ys)
829	
830	def counted_at(+n: U32, xs: List<&2, F32>, ys: List<&2, F32>) -> Series:
831	  Series{n, column(n, ys), At{column(n, xs)}}
832	
833	# The series of the values `ys` at `xs`, which must hold as many values,
834	# increasing.
835	def series_at(xs: List<&2, F32>, +ys: List<&2, F32>) -> Series:
836	  counted_at(length(ys, 0), xs, ys)
837	
838	def chosen(is_visibility: Bool, +g: Visibility, +dir: Direction, +k: Weight, +lo: F32, +hi: F32, +limit: U32, ys: Array<F32>, xs: Xs, +n: U32) -> Edges:
839	  match is_visibility:
840	    case True{}:
841	      visibility(g, dir, k, lo, hi, ys, xs, n)
842	    case False{}:
843	      penetrable(g, dir, k, lo, hi, limit, ys, xs, n)
844	
845	# The graph of `s` with the options `o`, with the edges ts2vg's `build`
846	# makes. The natural and horizontal graphs run on one thread, and the
847	# limited penetrable ones split their sweeps among every thread.
848	def build(o: Options, s: Series) -> Graph:
849	  Options{g, +dir, k, lo, hi, +limit} = o
850	  Series{+n, ys, xs} = s
851	  Graph{n, dir, chosen(U32.is_zero(limit), g, dir, k, lo, hi, limit, ys, xs, n)}
852	
853	# What to read off the edges
854	
855	# `f(acc, src, dst, w)` over every edge of `g`, in the order the rope
856	# holds them, from `acc`.
857	def fold.edges(~A: Type, ~f: A -> U32 -> U32 -> F32 -> A, es: Edges, acc: A) -> A:
858	  match es:
859	    case ENil{}:
860	      acc
861	    case ECon{s, t, w, next}:
862	      fold.edges(~A, ~f, next, f(acc, s, t, w))
863	    case EJoin{a, b}:
864	      fold.edges(~A, ~f, b, fold.edges(~A, ~f, a, acc))
865	
866	def fold(~A: Type, ~f: A -> U32 -> U32 -> F32 -> A, g: Graph, acc: A) -> A:
867	  Graph{n, dir, es} = g
868	  fold.edges(~A, ~f, es, acc)
869	
870	# The degree of each node of `g`, as ts2vg's `degrees`: the edges at the
871	# node, in either direction.
872	def degrees(g: Graph) -> Array<U32>:
873	  Graph{+n, dir, es} = g
874	  degrees.into(es, slots(n))
875	
876	# Each half edge (u, v), from u to v, goes into two passes of a counting
877	# sort: by v, then by u, so each row ends up in increasing order.
878	
879	def is_undirected(d: Direction) -> Bool:
880	  match d:
881	    case Undirected{}:
882	      True{}
883	    case _:
884	      False{}
885	
886	def bumped(r: Array<U32> & U32, +i: U32) -> Array<U32>:
887	  (a, +c) = r
888	  Array.set(U32, a, i, (c + 1 : U32))
889	
890	type Counts is Type:
891	  Counts{outs: Array<U32>, ins: Array<U32>}
892	
893	def count_half(present: Bool, +u: U32, +v: U32, c: Counts) -> Counts:
894	  match present:
895	    case True{}:
896	      Counts{outs, ins} = c
897	      Counts{bumped(Array.get(U32, outs, u), u), bumped(Array.get(U32, ins, v), v)}
898	    case False{}:
899	      c
900	
901	def count_halves(es: Edges, +both: Bool, c: Counts) -> Counts:
902	  match es:
903	    case ENil{}:
904	      c
905	    case ECon{+s, +t, w, next}:
906	      count_halves(next, both, count_half(both, t, s, count_half(True{}, s, t, c)))
907	    case EJoin{a, b}:
908	      count_halves(b, both, count_halves(a, both, c))
909	
910	# The counts of `n` nodes, in place, as where each node's run starts; slot
911	# `n` gets the total.
912	def starts(count: Nat, +i: U32, r: Array<U32> & U32, +acc: U32) -> Array<U32>:
913	  match count:
914	    case 0n:
915	      (a, c) = r
916	      Array.set(U32, a, i, acc)
917	    case 1n+k:
918	      (a, +c) = r
919	      +i2 = (i + 1 : U32)
920	      starts(k, i2, Array.get(U32, Array.set(U32, a, i, acc), i2), (acc + c : U32))
921	
922	# The half edges, by target: source, target and weight at each position.
923	type Halves is Type:
924	  Halves{cur: Array<U32>, us: Array<U32>, vs: Array<U32>, ws: Array<F32>}
925	
926	def put_half(+u: U32, +v: U32, +w: F32, r: Array<U32> & U32, us: Array<U32>, vs: Array<U32>, ws: Array<F32>) -> Halves:
927	  (cur, +p) = r
928	  Halves{Array.set(U32, cur, v, (p + 1 : U32)), Array.set(U32, us, p, u), Array.set(U32, vs, p, v), Array.set(F32, ws, p, w)}
929	
930	def place_half(present: Bool, +u: U32, +v: U32, +w: F32, h: Halves) -> Halves:
931	  match present:
932	    case True{}:
933	      Halves{cur, us, vs, ws} = h
934	      put_half(u, v, w, Array.get(U32, cur, v), us, vs, ws)
935	    case False{}:
936	      h
937	
938	def by_target(es: Edges, +both: Bool, h: Halves) -> Halves:
939	  match es:
940	    case ENil{}:
941	      h
942	    case ECon{+s, +t, +w, next}:
943	      by_target(next, both, place_half(both, t, s, w, place_half(True{}, s, t, w, h)))
944	    case EJoin{a, b}:
945	      by_target(b, both, by_target(a, both, h))
946	
947	type Rows is Type:
948	  Rows{cur: Array<U32>, targets: Array<U32>, weights: Array<F32>}
949	
950	def put_row(+u: U32, +v: U32, +w: F32, r: Array<U32> & U32, targets: Array<U32>, weights: Array<F32>) -> Rows:
951	  (cur, +p) = r
952	  Rows{Array.set(U32, cur, u, (p + 1 : U32)), Array.set(U32, targets, p, v), Array.set(F32, weights, p, w)}
953	
954	# The half edges from position `p` on, in target order, into their
955	# sources' rows.
956	def by_source(count: Nat, +p: U32, ur: Array<U32> & U32, vr: Array<U32> & U32, wr: Array<F32> & F32, out: Rows) -> Rows:
957	  match count:
958	    case 0n:
959	      (us, u) = ur
960	      (vs, v) = vr
961	      (ws, w) = wr
962	      out
963	    case 1n+k:
964	      (us, +u) = ur
965	      (vs, +v) = vr
966	      (ws, +w) = wr
967	      Rows{cur, targets, weights} = out
968	      +p2 = (p + 1 : U32)
969	      by_source(k, p2, Array.get(U32, us, p2), Array.get(U32, vs, p2), Array.get(F32, ws, p2), put_row(u, v, w, Array.get(U32, cur, u), targets, weights))
970	
971	def slots_for(+n: U32) -> Nat:
972	  Nat.add(1n, U32.log2(n))
973	
974	def rows_done(offsets: Array<U32>, r: Rows) -> Adjacency:
975	  Rows{cur, targets, weights} = r
976	  Adjacency{offsets, targets, weights}
977	
978	def rows_from(+h: U32, oc: Array<U32> & Array<U32>, hv: Halves) -> Adjacency:
979	  (offsets, cur) = oc
980	  Halves{tcur, us, vs, ws} = hv
981	  rows_done(offsets, by_source(U32.to_nat(h), 0, Array.get(U32, us, 0), Array.get(U32, vs, 0), Array.get(F32, ws, 0), Rows{cur, Array.new(U32, slots_for(h), 0), Array.new(F32, slots_for(h), 0.0)}))
982	
983	def with_starts(es: Edges, +both: Bool, +h: U32, outs: Array<U32>, ins: Array<U32>) -> Adjacency:
984	  rows_from(h, Array.clone(U32, outs), by_target(es, both, Halves{ins, Array.new(U32, slots_for(h), 0), Array.new(U32, slots_for(h), 0), Array.new(F32, slots_for(h), 0.0)}))
985	
986	def with_total(es: Edges, +both: Bool, r: Array<U32> & U32, ins: Array<U32>) -> Adjacency:
987	  (outs, +h) = r
988	  with_starts(es, both, h, outs, ins)
989	
990	def tallied(es: Edges, +both: Bool, +n: U32, c: Counts) -> Adjacency:
991	  Counts{outs, ins} = c
992	  with_total(es, both, Array.get(U32, starts(U32.to_nat(n), 0, Array.get(U32, outs, 0), 0), n), starts(U32.to_nat(n), 0, Array.get(U32, ins, 0), 0))
993	
994	# The rows of the nodes of `g`.
995	def adjacency(g: Graph) -> Adjacency:
996	  Graph{+n, +dir, +es} = g
997	  tallied(es, is_undirected(dir), n, count_halves(es, is_undirected(dir), Counts{Array.new(U32, slots_for(n), 0), Array.new(U32, slots_for(n), 0)}))