src/penetrable.bend
bend-visgraph@1.0.1.0 · 228 lines · 0 definitions
1# Limited penetrable visibility graphs, as ts2vg 1.2.4's _natural_penetrable
2# and _horizontal_penetrable build them, in F32: every point sweeps the
3# points right of it, and an edge may cross up to `limit` points. The sweeps
4# are independent, so the points split into parallel ranges of equal work.
5import Base
6import bend-visgraph@1.0.1.0/src/graph.bend as G
7
8# The points, in a list the parallel ranges share.
9type Pts is Data:
10 PNil{}
11 PCon{i: U32, x: F32, y: F32, next: Pts}
12
13# The points of `ys` and `xs` from index `i` down to 0, onto `acc`; `yr`
14# and `xr` hold y and x at `i`, and x is the index if `index` says so. The
15# last read, at index -1, wraps around.
16def points(count: Nat, +i: U32, yr: Array<F32> & F32, xr: Array<F32> & F32, +index: Bool, acc: Pts) -> Pts:
17 match count:
18 case 0n:
19 (ys, y) = yr
20 (xs, x) = xr
21 acc
22 case 1n+c:
23 (ys, +y) = yr
24 (xs, +x) = xr
25 +i2 = (i - 1 : U32)
26 points(c, i2, Array.get(F32, ys, i2), Array.get(F32, xs, i2), index, PCon{i, G.pick(index, U32.to_f32(i), x), y, acc})
27
28# The `limit + 1` largest slopes (or heights) a sweep has accepted, as ts2vg's
29# max_slopes; the sweep accepts a point above the smallest of them.
30type Tops is Data:
31 TNil{}
32 TCon{v: F32, next: Tops}
33
34# The position and value of the first smallest top.
35type Low is Data:
36 Low{k: Nat, v: F32}
37
38def tops_new(n: Nat) -> Tops:
39 match n:
40 case 0n:
41 TCon{G.neg_inf(), TNil{}}
42 case 1n+p:
43 TCon{G.neg_inf(), tops_new(p)}
44
45def pick_low(c: Bool, a: Low, b: Low) -> Low:
46 match c:
47 case True{}:
48 a
49 case False{}:
50 b
51
52def low_from(ts: Tops, +k: Nat, +best: Low) -> Low:
53 match ts:
54 case TNil{}:
55 best
56 case TCon{+v, next}:
57 Low{bk, +bv} = best
58 low_from(next, 1n+k, pick_low(F32.is_lt(v, bv), Low{k, v}, Low{bk, bv}))
59
60def lowest(ts: Tops) -> Low:
61 low_from(ts, 0n, Low{0n, G.inf()})
62
63def Low.v(+l: Low) -> F32:
64 Low{k, v} = l
65 v
66
67def set_at(k: Nat, ts: Tops, v: F32) -> Tops:
68 match k ts:
69 case 0n TCon{old, next}:
70 TCon{v, next}
71 case 1n+p TCon{old, next}:
72 TCon{old, set_at(p, next, v)}
73 case _ TNil{}:
74 TNil{}
75
76# How many tops a natural edge with slope `s` does not clear, as ts2vg's
77# num_penetrations weight.
78def blocked_natural(ts: Tops, +s: F32, +t: F32, +n: F32) -> F32:
79 match ts:
80 case TNil{}:
81 n
82 case TCon{v, next}:
83 blocked_natural(next, s, t, G.pick(G.greater(s, v, t), n, F32.add(n, 1.0)))
84
85def blocked_horizontal(ts: Tops, +ya: F32, +yb: F32, +n: F32) -> F32:
86 match ts:
87 case TNil{}:
88 n
89 case TCon{+v, next}:
90 blocked_horizontal(next, ya, yb, G.pick(F32.is_le(ya, v) || F32.is_le(yb, v), F32.add(n, 1.0), n))
91
92# ts2vg's add_edge with a known weight: NaN means the weight function
93# computes it.
94def add_known(+k: G.Weight, +lo: F32, +hi: F32, i1: U32, i2: U32, +x1: F32, +x2: F32, +y1: F32, +y2: F32, +s: F32, +known: F32, es: G.Edges) -> G.Edges:
95 +w = G.pick(G.is_nan(known), G.weight(k, x1, x2, y1, y2, s), known)
96 G.put(F32.is_le(w, lo) || F32.is_ge(w, hi), i1, i2, w, es)
97
98# The edge from a to b, or from b to a for top_to_bottom when b is higher.
99def add_oriented(down: Bool, +k: G.Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, xa: F32, xb: F32, ya: F32, yb: F32, s: F32, known: F32, es: G.Edges) -> G.Edges:
100 match down:
101 case True{}:
102 add_known(k, lo, hi, ib, ia, xb, xa, yb, ya, s, known, es)
103 case False{}:
104 add_known(k, lo, hi, ia, ib, xa, xb, ya, yb, s, known, es)
105
106def is_down(dir: G.Direction, +ya: F32, +yb: F32) -> Bool:
107 match dir:
108 case G.TopToBottom{}:
109 F32.is_gt(yb, ya)
110 case _:
111 False{}
112
113def is_count(k: G.Weight) -> Bool:
114 match k:
115 case G.NumPenetrations{}:
116 True{}
117 case _:
118 False{}
119
120def natural_edge(+dir: G.Direction, +k: G.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: G.Edges) -> G.Edges:
121 +known = G.pick(is_count(k), blocked_natural(ts, s, G.tol(scale_a, xb, yb), 0.0), G.nan())
122 add_oriented(is_down(dir, ya, yb), k, lo, hi, ia, ib, xa, xb, ya, yb, s, known, es)
123
124def horizontal_edge(+dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, ia: U32, ib: U32, +xa: F32, xb: F32, +ya: F32, +yb: F32, ts: Tops, es: G.Edges) -> G.Edges:
125 +known = G.pick(is_count(k), blocked_horizontal(ts, ya, yb, 0.0), G.nan())
126 add_oriented(is_down(dir, ya, yb), k, lo, hi, ia, ib, xa, xb, ya, yb, G.nan(), known, es)
127
128# The sweep of point a over the points right of it: ts2vg's inner loop. A
129# point b the sweep sees replaces the smallest top and adds its edge at the
130# next step, when `up` says so, so a step after an unseen point calls no
131# def that recurses.
132def natural_sweep(pts: Pts, up: Bool, +dir: G.Direction, +k: G.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: G.Edges) -> G.Edges:
133 match pts up:
134 case PNil{} False{}:
135 es
136 case PNil{} True{}:
137 natural_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, sb, scale_a, ts, es)
138 case PCon{+ic, +xc, +yc, next} False{}:
139 Low{+lk, +lv} = low
140 +s = F32.div(F32.sub(yc, ya), F32.sub(xc, xa))
141 natural_sweep(next, G.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)
142 case PCon{+ic, +xc, +yc, next} True{}:
143 Low{+lk, lv} = low
144 +after = set_at(lk, ts, sb)
145 +now = lowest(after)
146 +s = F32.div(F32.sub(yc, ya), F32.sub(xc, xa))
147 natural_sweep(next, G.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))
148
149# The same for the horizontal sweep, which stops once the smallest top is
150# above a: no later point can be seen.
151def horizontal_sweep(pts: Pts, up: Bool, stop: Bool, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +ia: U32, +xa: F32, +ya: F32, +ts: Tops, low: Low, +ib: U32, +xb: F32, +yb: F32, es: G.Edges) -> G.Edges:
152 match pts up stop:
153 case _ False{} True{}:
154 es
155 case PNil{} False{} False{}:
156 es
157 case PNil{} True{} _:
158 horizontal_edge(dir, k, lo, hi, ia, ib, xa, xb, ya, yb, ts, es)
159 case PCon{+ic, xc, +yc, next} False{} False{}:
160 Low{+lk, +lv} = low
161 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)
162 case PCon{+ic, xc, +yc, next} True{} _:
163 Low{+lk, lv} = low
164 +after = set_at(lk, ts, yb)
165 +now = lowest(after)
166 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))
167
168def sweep(g: G.Graph, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +limit: Nat, ia: U32, +xa: F32, +ya: F32, pts: Pts, es: G.Edges) -> G.Edges:
169 match g:
170 case G.Natural{}:
171 natural_sweep(pts, False{}, dir, k, lo, hi, ia, xa, ya, G.scale(xa, ya), tops_new(limit), Low{0n, G.neg_inf()}, 0, 0.0, 0.0, 0.0, es)
172 case G.Horizontal{}:
173 horizontal_sweep(pts, False{}, False{}, dir, k, lo, hi, ia, xa, ya, tops_new(limit), Low{0n, G.neg_inf()}, 0, 0.0, 0.0, es)
174
175# The sweeps of the first `count` points of `pts`, one rope each: ts2vg's
176# outer loop.
177def sweeps(pts: Pts, count: Nat, +g: G.Graph, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +limit: Nat) -> G.Edges:
178 match pts count:
179 case PCon{i, x, y, +next} 1n+c:
180 G.EJoin{sweep(g, dir, k, lo, hi, limit, i, x, y, next, G.ENil{}), sweeps(next, c, g, dir, k, lo, hi, limit)}
181 case _ _:
182 G.ENil{}
183
184def drop(n: Nat, pts: Pts) -> Pts:
185 match n pts:
186 case 1n+p PCon{i, x, y, next}:
187 drop(p, next)
188 case _ _:
189 pts
190
191# How many of the first `count` of `len` points do half the sweeps' work. A
192# natural sweep walks every point right of its start, so the first c starts
193# walk c (2 len - 1 - c) / 2 points; a horizontal sweep stops early, so
194# halving the starts is close enough.
195def first_half(g: G.Graph, +len: U32, +count: U32) -> U32:
196 match g:
197 case G.Natural{}:
198 +b = F32.sub(F32.mul(2.0, U32.to_f32(len)), 1.0)
199 +c = U32.to_f32(count)
200 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)))
201 case G.Horizontal{}:
202 U32.div(count, 2)
203
204# The sweeps of the first `count` points of `pts`, as two halves of equal
205# work in parallel, `depth` times over. `len` and `cu` are `pts`'s length
206# and `count` as U32; they only place the cut.
207def split(depth: Nat, +g: G.Graph, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +limit: Nat, +pts: Pts, +len: U32, +cu: U32, +count: Nat) -> G.Edges:
208 match depth:
209 case 0n:
210 sweeps(pts, count, g, dir, k, lo, hi, limit)
211 case 1n++d:
212 +c1u = first_half(g, len, cu)
213 +c1 = Nat.min(count, U32.to_nat(c1u))
214 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))
215 G.EJoin{a, b}
216
217def built(+index: Bool, +g: G.Graph, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +limit: U32, ys: Array<F32>, xs: Array<F32>, +n: U32) -> G.Edges:
218 split(7n, g, dir, k, lo, hi, U32.to_nat(limit), 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))
219
220# The limited penetrable graph of the `n` values `ys` at `xs`, crossing at
221# most `limit` points. The x of `G.Index{}` is the index as an F32, exact
222# for up to 2^24 values; the x of `G.At{xs}` must increase.
223def build(+g: G.Graph, +dir: G.Direction, +k: G.Weight, +lo: F32, +hi: F32, +limit: U32, ys: Array<F32>, xs: G.Xs, +n: U32) -> G.Edges:
224 match xs:
225 case G.Index{}:
226 built(True{}, g, dir, k, lo, hi, limit, ys, G.no_xs(), n)
227 case G.At{a}:
228 built(False{}, g, dir, k, lo, hi, limit, ys, a, n)