Merge pull request #333 from tdamsma/perf/basis-eval-stack

perf: point evaluation without a heap basis array (outputs bit-identical)
This commit is contained in:
bherbruck
2026-10-21 14:38:45 +00:00
committed by GitHub
4 changed files with 256 additions and 28 deletions
+4 -4
View File
@@ -286,15 +286,15 @@ impl Curve {
return Err(SolidError::nlib(OP, NlibCode::Parameter));
}
// Compute non-vanishing B-splines.
let basis = knt.basis_eval(self.degree(), u, side)?;
let spn = sx(basis.span, OP)?;
// Compute non-vanishing B-splines (on the stack, as upstream's `N`).
let mut n = [0.0_f64; MAX_DEGREE + 2];
let spn = sx(knt.basis_eval_into(self.degree(), u, side, &mut n)?, OP)?;
// Compute the point on the curve.
let mut cw = CPoint::new(0.0, 0.0, 0.0, 0.0);
let j = spn - p;
for i in j..=spn {
let alpha = get(&basis.values, i - j, OP)?;
let alpha = get(&n, i - j, OP)?;
let q = get(pw, i, OP)?;
cw.x += alpha * q.x;
cw.y += alpha * q.y;
+66 -17
View File
@@ -18,6 +18,7 @@
//! | [`KnotVector::find_span_and_mult`] | `N_BasisFindSpanAndMult` |
//! | [`KnotVector::span_count`] | `N_BasisGetSpanCount` |
//! | [`KnotVector::basis_eval`] | `N_BasisEval` |
//! | [`KnotVector::basis_eval_into`] | `N_BasisEval` (caller's `N` array) |
//! | [`KnotVector::basis_i_eval`] | `N_BasisIEval` |
//! | [`KnotVector::basis_i_derivs`] | `N_BasisIDerivs` |
//! | [`Side`] | `NL_LEFT` / `NL_RIGHT` flags |
@@ -112,6 +113,14 @@ fn out_of_range(operation: &'static str) -> SolidError {
}
}
/// `p > NL_MAXDEG + 1`: upstream's `N_BasisEval` scratch arrays would overflow.
fn degree_too_high(operation: &'static str) -> SolidError {
SolidError::InvalidInput {
operation,
reason: "degree above NL_MAXDEG + 1 (upstream scratch array size)",
}
}
/// `a[i]`, or `InvalidInput` when out of range.
#[inline]
fn rd(a: &[f64], i: usize, op: &'static str) -> Result<f64> {
@@ -418,6 +427,20 @@ impl<'a> KnotVector<'a> {
/// (more than `p + 1` equal knots); `InvalidInput` where upstream would index out
/// of bounds (fewer than `2p + 2` knots, `p > MAX_DEGREE + 1`).
pub fn basis_eval(&self, p: usize, u: f64, side: Side) -> Result<Basis> {
const OP: &str = OP_BASIS_EVAL;
if p > MAX_DEGREE + 1 {
return Err(degree_too_high(OP));
}
let mut n = zeros(Some(p + 1), OP)?;
let span = self.basis_eval_into(p, u, side, &mut n)?;
Ok(Basis { span, values: n })
}
/// [`KnotVector::basis_eval`] into a caller buffer: zeroes `n[0..=p]`, writes the
/// values there and returns the span (same values and errors; `InvalidInput` when
/// `n` has fewer than `p + 1` entries). Lets hot evaluators keep `N` on the stack
/// (`[f64; MAX_DEGREE + 2]` fits every valid degree), as upstream does.
pub fn basis_eval_into(&self, p: usize, u: f64, side: Side, n: &mut [f64]) -> Result<usize> {
const OP: &str = OP_BASIS_EVAL;
// Divide by zero? The denominator is the sum of the two values multiplied
// by it, so the quotients never exceed 1; only an exactly-zero delta (more
@@ -426,25 +449,19 @@ impl<'a> KnotVector<'a> {
let uk = self.knots;
let m = self.last_index();
if p > MAX_DEGREE + 1 {
return Err(SolidError::InvalidInput {
operation: OP,
reason: "degree above NL_MAXDEG + 1 (upstream scratch array size)",
});
return Err(degree_too_high(OP));
}
let mut n = zeros(Some(p + 1), OP)?;
n.get_mut(..=p).ok_or_else(|| out_of_range(OP))?.fill(0.0);
// Special cases for end values.
if u == rd(uk, p, OP)? {
wr(&mut n, 0, 1.0, OP)?;
return Ok(Basis { span: p, values: n });
wr(n, 0, 1.0, OP)?;
return Ok(p);
}
let end = sub(m, p, OP)?;
if u == rd(uk, end, OP)? {
wr(&mut n, p, 1.0, OP)?;
return Ok(Basis {
span: sub(end, 1, OP)?,
values: n,
});
wr(n, p, 1.0, OP)?;
return sub(end, 1, OP);
}
// Find the knot span u is in.
@@ -453,7 +470,7 @@ impl<'a> KnotVector<'a> {
// Compute the non-vanishing B-splines.
let mut left = [0.0_f64; MAX_DEGREE + 2];
let mut right = [0.0_f64; MAX_DEGREE + 2];
wr(&mut n, 0, 1.0, OP)?;
wr(n, 0, 1.0, OP)?;
for k in 1..=p {
wr(&mut left, k, u - rd(uk, sub(i + 1, k, OP)?, OP)?, OP)?;
wr(&mut right, k, rd(uk, i + k, OP)? - u, OP)?;
@@ -468,15 +485,15 @@ impl<'a> KnotVector<'a> {
code: NlibCode::Numerical,
});
}
let temp = rd(&n, l, OP)?;
let temp = rd(n, l, OP)?;
// Keep the divisions inside the parentheses: dividing N[l] by the
// knot difference can be huge.
wr(&mut n, l, saved + temp * (r / delta), OP)?;
wr(n, l, saved + temp * (r / delta), OP)?;
saved = temp * (lf / delta);
}
wr(&mut n, k, saved, OP)?;
wr(n, k, saved, OP)?;
}
Ok(Basis { span: i, values: n })
Ok(i)
}
/// The single basis function `N[i, p](u)` (upstream `N_BasisIEval`).
@@ -842,6 +859,38 @@ mod tests {
assert_eq!(start.values, vec![1.0, 0.0, 0.0]);
}
#[test]
fn basis_eval_into_matches_basis_eval_bitwise() {
let k = kv(&EX23);
for s in -5..=55 {
let u = f64::from(s) * 0.1;
for side in [Side::Left, Side::Right] {
for p in 0..=2 {
// a dirty buffer: every used entry must be overwritten
let mut n = [7.0_f64; MAX_DEGREE + 2];
let r = k.basis_eval_into(p, u, side, &mut n);
match k.basis_eval(p, u, side) {
Ok(b) => {
assert_eq!(r.unwrap(), b.span);
let bits: Vec<u64> = b.values.iter().map(|x| x.to_bits()).collect();
let got: Vec<u64> = n[..=p].iter().map(|x| x.to_bits()).collect();
assert_eq!(got, bits, "u={u} p={p}");
}
Err(e) => assert_eq!(r.unwrap_err(), e),
}
}
}
}
let mut short = [0.0; 2];
assert!(k.basis_eval_into(2, 0.5, Side::Left, &mut short).is_err());
let mut n = [0.0; MAX_DEGREE + 3];
assert_eq!(
k.basis_eval_into(MAX_DEGREE + 2, 0.5, Side::Left, &mut n)
.unwrap_err(),
k.basis_eval(MAX_DEGREE + 2, 0.5, Side::Left).unwrap_err()
);
}
#[test]
fn basis_eval_partition_of_unity() {
let k = kv(&EX23);
+9 -7
View File
@@ -727,15 +727,17 @@ impl Surface {
let (pu, qu) = self.degrees();
let (knu, knv) = self.knot_vectors()?;
// Compute non-vanishing B-splines.
let bu = knu.basis_eval(pu, u, ufl)?;
let bv = knv.basis_eval(qu, v, vfl)?;
// Compute non-vanishing B-splines (on the stack, as upstream's `Nu`/`Nv`).
let mut nu = [0.0_f64; MAX_DEGREE + 2];
let mut nv = [0.0_f64; MAX_DEGREE + 2];
let u_span = knu.basis_eval_into(pu, u, ufl, &mut nu)?;
let v_span = knv.basis_eval_into(qu, v, vfl, &mut nv)?;
// Compute the point on the surface.
let (p, q) = (sx(pu, OP)?, sx(qu, OP)?);
let mut sw = CPoint::ZERO;
let tmp1 = sx(bu.span, OP)? - p;
let tmp2 = sx(bv.span, OP)? - q;
let tmp1 = sx(u_span, OP)? - p;
let tmp2 = sx(v_span, OP)? - q;
for i in 0..=p {
let mut tw = CPoint::ZERO;
let row = self.row(ux(tmp1 + i, OP)?)?;
@@ -743,12 +745,12 @@ impl Surface {
blend_spec(
pnoz,
pnow,
get(&bv.values, j, OP)?,
get(&nv, j, OP)?,
get(row, tmp2 + j, OP)?,
&mut tw,
);
}
blend_spec(pnoz, pnow, get(&bu.values, i, OP)?, tw, &mut sw);
blend_spec(pnoz, pnow, get(&nu, i, OP)?, tw, &mut sw);
}
Ok(to_euclid_spec(pnoz, pnow, sw))
}
+177
View File
@@ -0,0 +1,177 @@
#!/usr/bin/env python3
"""Top-function tables from a samply profile (Firefox processed format), no GUI needed.
usage: samply_top.py PROFILE.json.gz [--top N] [--focus SUBSTR] [--thread NAME] [--callees SUBSTR] [--alloc]
Record with `samply record --save-only --unstable-presymbolicate -o P.json.gz -- <cmd>`.
The main binary's frames are expanded to inline level with gimli `addr2line -i`
(`cargo install addr2line --features bin`); other libraries use the `.syms.json`
sidecar (function level). Prints self (innermost frame) and inclusive (once per
sample) percentages.
"""
import bisect
import collections
import gzip
import json
import os
import re
import subprocess
import sys
def short(name):
name = re.sub(r"::h[0-9a-f]{16}$", "", name)
name = re.sub(r"<impl [^>]*>", "<impl>", name)
return name if len(name) < 150 else name[:147] + "..."
def arg(flag, default):
return sys.argv[sys.argv.index(flag) + 1] if flag in sys.argv else default
def main():
path = sys.argv[1]
top = int(arg("--top", "40"))
focus = arg("--focus", None)
want_thread = arg("--thread", None)
prof = json.load(gzip.open(path) if path.endswith(".gz") else open(path))
syms_path = re.sub(r"\.gz$", "", path) + ".syms.json"
if not os.path.exists(syms_path):
syms_path = path.replace(".json.gz", ".json.syms.json")
syms = json.load(open(syms_path)) if os.path.exists(syms_path) else {"string_table": [], "data": []}
st = syms["string_table"]
tables = {}
for e in syms["data"]:
rows = sorted((r["rva"], r["size"], r["symbol"]) for r in e["symbol_table"])
tables[e["debug_name"]] = ([r[0] for r in rows], rows)
libs = prof["libs"]
main_lib = None
# samples: list of frame lists (outermost first) as (libname, rva, is_leaf)
stacks = collections.Counter()
for t in prof["threads"]:
if want_thread and t["name"] != want_thread:
continue
ft, stt, res = t["frameTable"], t["stackTable"], t["resourceTable"]
funcs = t["funcTable"]
memo = {}
def frames(s):
out = []
while s is not None:
f = stt["frame"][s]
func = ft["func"][f]
r = funcs["resource"][func]
lib = libs[res["lib"][r]]["debugName"] if r is not None and r >= 0 and res["lib"][r] is not None else "?"
out.append((lib, ft["address"][f]))
s = stt["prefix"][s]
return tuple(reversed(out))
smp = t["samples"]
weights = smp.get("weight") or [1] * smp["length"]
for s, w in zip(smp["stack"], weights):
if s is None:
continue
if s not in memo:
memo[s] = frames(s)
stacks[memo[s]] += w or 1
total = sum(stacks.values())
exe = next((l for l in libs if l["debugName"] == "difftest"), None)
addrs = set()
for fr, _ in stacks.items():
for i, (lib, a) in enumerate(fr):
if exe and lib == exe["debugName"] and a >= 0:
addrs.add(a if i == len(fr) - 1 else a - 1)
sym = {}
if exe and addrs:
alist = sorted(addrs)
proc = subprocess.run(
[os.path.expanduser("~/.cargo/bin/addr2line"), "-f", "-C", "-i", "-a", "-e", exe["path"]] + [hex(a) for a in alist],
capture_output=True, text=True)
out = proc.stdout.splitlines()
cur = None
i = 0
while i < len(out):
l = out[i]
if l.startswith("0x"):
cur = int(l, 16)
sym[cur] = []
i += 1
continue
sym[cur].append(short(l))
i += 2
def name_of(lib, a, leaf):
if exe and lib == exe["debugName"]:
k = a if leaf else a - 1
chain = sym.get(k)
if chain:
return list(reversed(chain)) # outermost first
tb = tables.get(lib)
if tb:
j = bisect.bisect_right(tb[0], a) - 1
if j >= 0:
rva, size, si = tb[1][j]
if a < rva + max(size, 1) and si < len(st):
return [short(st[si]) + " [" + lib + "]"]
return ["[" + lib + "]"]
selfc = collections.Counter()
incl = collections.Counter()
split = collections.Counter()
callees = collections.Counter()
allocs = collections.Counter()
alloc_mode = "--alloc" in sys.argv
callees_of = arg("--callees", None)
n = 0
for fr, w in stacks.items():
names = []
for i, (lib, a) in enumerate(fr):
names.extend(name_of(lib, a, i == len(fr) - 1))
inner = names[-1] if names else "?"
cxx = "(" in inner.split("<")[0] or "libstdc++" in inner
if any("::run_cpp" in x for x in names):
split["run_cpp"] += w
elif cxx or any(x.startswith("solidcore_ffi::") for x in names):
split["C++ outside run_cpp (FFI calls from run_rust, or truncated C++ stacks)"] += w
elif any(x.startswith("solidcore::") or x.startswith("<solidcore::") for x in names):
split["Rust kernel (solidcore::)"] += w
elif any("difftest::" in x for x in names):
split["harness (difftest)"] += w
else:
split["libc / unattributed"] += w
if focus and not any(focus in x for x in names):
continue
n += w
if callees_of:
idx = [i for i, x in enumerate(names) if callees_of in x]
if idx:
j = idx[-1] + 1
callees[names[j] if j < len(names) else "(self)"] += w
selfc[names[-1] if names else "?"] += w
if alloc_mode and names and re.search(r"malloc|free|realloc|calloc|memset|memcpy|memmove|finish_grow|deallocate|try_allocate_in|__rust_alloc|__rust_dealloc|__rust_realloc|alloc_shim|fun_[0-9a-f]+ \[libc", names[-1]):
owner = next((x for x in reversed(names) if x.startswith(("solidcore::", "<solidcore::", "difftest::", "<difftest::")) or re.match(r"^[A-Za-z_]+\w*\(", x)), "?")
allocs[owner] += w
for x in set(names):
incl[x] += w
print(f"samples: {total} total, {n} selected" + (f" (focus {focus!r})" if focus else ""))
print("\n## split (all samples)\n| % | samples | where |\n|---:|---:|---|")
for k, v in sorted(split.items(), key=lambda kv: -kv[1]):
print(f"| {100*v/max(total,1):.1f} | {v} | {k} |")
if callees_of:
print(f"\n## callees of {callees_of!r} (next frame below its innermost occurrence)\n| % | samples | callee |\n|---:|---:|---|")
tot = sum(callees.values())
for k, v in callees.most_common(top):
print(f"| {100*v/max(tot,1):.1f} | {v} | `{k}` |")
if alloc_mode:
print(f"\n## allocator / memset samples by nearest Rust (or C++) owner\n| % of selected | samples | owner |\n|---:|---:|---|")
for k, v in allocs.most_common(top):
print(f"| {100*v/max(n,1):.1f} | {v} | `{k}` |")
print(f"\n## self (top {top})\n| % | samples | function |\n|---:|---:|---|")
for k, v in selfc.most_common(top):
print(f"| {100*v/max(n,1):.1f} | {v} | `{k}` |")
print(f"\n## inclusive (top {top})\n| % | samples | function |\n|---:|---:|---|")
for k, v in incl.most_common(top):
print(f"| {100*v/max(n,1):.1f} | {v} | `{k}` |")
if __name__ == "__main__":
main()