visgraph.bend

bend-visgraph@2.0.0.0 · 997 lines · 134 definitions

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.
11import Base
12
13type Visibility is Data:
14  Natural{}
15  Horizontal{}
16
17type Direction is Data:
18  Undirected{}
19  LeftToRight{}
20  TopToBottom{}
21
22type 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.
38type 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.
44type 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.
53type 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.
58type 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.
62type 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.
70type 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
77def nan() -> F32:
78  F32.div(0.0, 0.0)
79
80def inf() -> F32:
81  F32.div(1.0, 0.0)
82
83def neg_inf() -> F32:
84  F32.div(F32.neg(1.0), 0.0)
85
86def is_nan(+a: F32) -> Bool:
87  F32.is_ne(a, a)
88
89def pick(c: Bool, a: F32, b: F32) -> F32:
90  match c:
91    case True{}:
92      a
93    case False{}:
94      b
95
96def 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
106def 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
109def 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
112def 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
142def 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
149def 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).
154def 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).
163def 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
170def 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
179def abs_tol() -> F32:
180  {0.00000000000001 : F32}
181
182def rel_tol() -> F32:
183  {0.00000000000001 : F32}
184
185# max(|xa|, |ya|), which a sweep computes once.
186def 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.
192def 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
195def 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.
201def 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.
217type 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.
226def 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.
230type Popped is Type:
231  Popped{first: U32, below: F32, st: Stack, rights: Array<U32>, visits: F32}
232
233def 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
244def 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.
253type Cols is Type:
254  Cols{ys: Array<F32>, lefts: Array<U32>, rights: Array<U32>}
255
256type 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
261def 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
272def 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.
284type 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.
290type Block is Data:
291  Block{start: U32, opening_nans: U32, unbounded: Unbounded, st: Stack, visits: F32}
292
293type Local is Type:
294  Local{b: Block, c: Cols}
295
296def 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
322def 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
331type Merged is Type:
332  Merged{st: Stack, lefts: Array<U32>, rights: Array<U32>, visits: F32}
333
334def 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
341def 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
348def 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
355def 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
360def 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
367def 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
375def 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
382def 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.
389def 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.
397type Walked is Type:
398  Walked{yr: Array<F32> & F32, xr: Array<F32> & F32, es: Edges}
399
400def at_index(xs: Array<F32>, +j: U32) -> Array<F32> & F32:
401  (xs, U32.to_f32(j))
402
403def 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
417def 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
431def 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
444def 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`.
458def 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
469def 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
488type Out is Type:
489  Out{cols: Cols, xs: Array<F32>, es: Edges}
490
491def 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
510def edges_of(o: Out) -> Edges:
511  Out{cols, xs, es} = o
512  es
513
514def 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
518def 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
522def 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.
528def 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.
532def no_xs() -> Array<F32>:
533  Array.new(F32, 0n, 0.0)
534
535def 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.
544def 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
554def added(r: Array<U32> & U32) -> Array<U32>:
555  (d, old) = r
556  d
557
558def 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.
563def 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.
578type 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.
585def 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.
599type Tops is Data:
600  TNil{}
601  TCon{v: F32, next: Tops}
602
603# The position and value of the first smallest top.
604type Low is Data:
605  Low{k: Nat, v: F32}
606
607def 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
614def pick_low(c: Bool, a: Low, b: Low) -> Low:
615  match c:
616    case True{}:
617      a
618    case False{}:
619      b
620
621def 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
629def lowest(ts: Tops) -> Low:
630  low_from(ts, 0n, Low{0n, inf()})
631
632def Low.v(+l: Low) -> F32:
633  Low{k, v} = l
634  v
635
636def 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.
647def 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
654def 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.
663def 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.
668def 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
675def 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
682def is_count(k: Weight) -> Bool:
683  match k:
684    case NumPenetrations{}:
685      True{}
686    case _:
687      False{}
688
689def 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
693def 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.
701def 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.
720def 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
737def 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.
746def 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
753def 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.
764def 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.
776def 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
786def 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.
792def 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.
803def options(v: Visibility) -> Options:
804  Options{v, Undirected{}, Unweighted{}, neg_inf(), inf(), 0}
805
806def 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
813def 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
820def 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
823def 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.
827def series(+ys: List<&2, F32>) -> Series:
828  counted(length(ys, 0), ys)
829
830def 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.
835def series_at(xs: List<&2, F32>, +ys: List<&2, F32>) -> Series:
836  counted_at(length(ys, 0), xs, ys)
837
838def 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.
848def 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`.
857def 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
866def 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.
872def 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
879def is_undirected(d: Direction) -> Bool:
880  match d:
881    case Undirected{}:
882      True{}
883    case _:
884      False{}
885
886def bumped(r: Array<U32> & U32, +i: U32) -> Array<U32>:
887  (a, +c) = r
888  Array.set(U32, a, i, (c + 1 : U32))
889
890type Counts is Type:
891  Counts{outs: Array<U32>, ins: Array<U32>}
892
893def 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
901def 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.
912def 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.
923type Halves is Type:
924  Halves{cur: Array<U32>, us: Array<U32>, vs: Array<U32>, ws: Array<F32>}
925
926def 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
930def 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
938def 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
947type Rows is Type:
948  Rows{cur: Array<U32>, targets: Array<U32>, weights: Array<F32>}
949
950def 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.
956def 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
971def slots_for(+n: U32) -> Nat:
972  Nat.add(1n, U32.log2(n))
973
974def rows_done(offsets: Array<U32>, r: Rows) -> Adjacency:
975  Rows{cur, targets, weights} = r
976  Adjacency{offsets, targets, weights}
977
978def 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
983def 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
986def 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
990def 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`.
995def 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)}))