crates/mimir-graph/src/ops.rs
//! The graph queries the lenses run (plan §6). `A` is the boolean adjacency
//! (`A[i][j]` = edge i → j, e.g. "i calls / depends on j").
//!
//! | op | algebra |
//! |---|---|
//! | [`reach`] | BFS, `q⟨¬v⟩ = q ANY.PAIR A` (vxm, push from the frontier) |
//! | [`reverse_reach`] | same on `Aᵀ` without materialising it: `q⟨¬v⟩ = A ANY.PAIR q` (mxv) |
//! | [`fan_out`] / [`fan_in`] | row reduce PLUS / `ones PLUS.PAIR A` |
//! | [`min_hops`] | level-synchronous MIN_PLUS (edge = 1 hop) from many seeds |
//! | [`max_on_path`] | fixpoint `s = max(s, A MAX.SECOND s)`: worst value reachable |
//! | [`bottleneck`] | (max, min) over NODE weights from seeds: `c = (f MAX.FIRST A) min w`, frontier fixpoint |
//! | [`mask_by`] / [`mask_vector`] | scope to a node set (diag(s)·A·diag(s), `w⟨s⟩ = v`) |
//! | [`similarity`] | `A·Aᵀ` with PLUS.PAIR (= PLUS.TIMES on a 0/1 pattern) |
//! | [`components`] | weakly connected components: MIN.SECOND label propagation + pointer jumping |
//! | [`diff`] | `added = Ay⟨¬Ax⟩`, `removed = Ax⟨¬Ay⟩` |
//! | [`mask_rows`] | `diag(set)·A`: only the out-edges of nodes in the set |
//! | [`union`] | `C = A ∨ B` (eWiseAdd LOR): several edge kinds walked as one graph |
//! | [`propagate`] | `y = A·x` over a WEIGHTED adjacency, PLUS.TIMES (FP64): one step of a random-projection embedding (FastRP) |
use std::ptr::null_mut;
use graphblas_sys as sys;
use crate::error::{Error, Result, grb};
use crate::matrix::Matrix;
use crate::vector::Vector;
fn square(a: &Matrix<bool>) -> Result<u64> {
let (r, c) = (a.nrows()?, a.ncols()?);
if r != c {
return Err(Error::invalid(format!(
"adjacency must be square, is {r}×{c}"
)));
}
Ok(r)
}
fn same_size<T: crate::Scalar>(v: &Vector<T>, n: u64, what: &str) -> Result<()> {
let s = v.size()?;
if s != n {
return Err(Error::invalid(format!(
"{what} has size {s}, graph has {n} nodes"
)));
}
Ok(())
}
/// Shared BFS loop. `forward`: push along out-edges (vxm); else along in-edges (mxv).
fn bfs(
seeds: &Vector<bool>,
a: &Matrix<bool>,
max_hops: Option<u32>,
forward: bool,
) -> Result<Vector<bool>> {
let n = square(a)?;
same_size(seeds, n, "seeds")?;
let visited = seeds.dup()?;
let frontier = seeds.dup()?;
let mut hops = 0u32;
loop {
if max_hops.is_some_and(|h| hops >= h) || frontier.nvals()? == 0 {
break;
}
// frontier⟨¬visited, struct, replace⟩ = next step from frontier
if forward {
grb!(GrB_vxm(
frontier.raw,
visited.raw,
null_mut(),
sys::GxB_ANY_PAIR_BOOL,
frontier.raw,
a.raw,
sys::GrB_DESC_RSC
))?;
} else {
grb!(GrB_mxv(
frontier.raw,
visited.raw,
null_mut(),
sys::GxB_ANY_PAIR_BOOL,
a.raw,
frontier.raw,
sys::GrB_DESC_RSC
))?;
}
// visited⟨frontier, struct⟩ = true
grb!(GrB_Vector_assign_BOOL(
visited.raw,
frontier.raw,
null_mut(),
true,
sys::GrB_ALL,
n,
sys::GrB_DESC_S
))?;
hops += 1;
}
Ok(visited)
}
/// Every node reachable from any seed in at most `max_hops` steps (`None` =
/// unbounded), seeds included. Values are `true`; absent = unreachable.
pub fn reach(
seeds: &Vector<bool>,
a: &Matrix<bool>,
max_hops: Option<u32>,
) -> Result<Vector<bool>> {
bfs(seeds, a, max_hops, true)
}
/// Every node that reaches a seed (the dependents walk / blast radius),
/// i.e. [`reach`] on `Aᵀ`, computed with mxv on `A` (no transpose is built).
/// For many walks on one graph, `reach(seeds, &a.transpose()?, ..)` with the
/// transpose kept around is faster — see the bench.
pub fn reverse_reach(
seeds: &Vector<bool>,
a: &Matrix<bool>,
max_hops: Option<u32>,
) -> Result<Vector<bool>> {
bfs(seeds, a, max_hops, false)
}
/// Out-degree per node (row reduce). Nodes with no out-edge are absent.
pub fn fan_out(a: &Matrix<bool>) -> Result<Vector<u64>> {
let w = Vector::<u64>::new(a.nrows()?)?;
grb!(GrB_Matrix_reduce_Monoid(
w.raw,
null_mut(),
null_mut(),
sys::GrB_PLUS_MONOID_UINT64,
a.raw,
null_mut()
))?;
Ok(w)
}
/// In-degree per node: `ones PLUS.PAIR A` (vxm). Measured ~4× faster than a
/// column reduce (`GrB_reduce` with a transposed input) on GraphBLAS 7.4.
/// Nodes with no in-edge are absent.
pub fn fan_in(a: &Matrix<bool>) -> Result<Vector<u64>> {
let n = a.nrows()?;
let ones = Vector::<u64>::new(n)?;
grb!(GrB_Vector_assign_UINT64(
ones.raw,
null_mut(),
null_mut(),
1,
sys::GrB_ALL,
n,
null_mut()
))?;
let w = Vector::<u64>::new(a.ncols()?)?;
grb!(GrB_vxm(
w.raw,
null_mut(),
null_mut(),
sys::GxB_PLUS_PAIR_UINT64,
ones.raw,
a.raw,
null_mut()
))?;
Ok(w)
}
/// Fewest hops from the nearest seed to every reachable node (seeds = 0),
/// bounded by `max_hops`. Level-synchronous MIN_PLUS: the frontier holds the
/// nodes first reached at hop k (all valued k), `t⟨¬dist⟩ = frontier MIN.PLUS A`
/// gives k+1 to the nodes not yet reached. With unit edges the first visit is
/// the minimum, so no relaxation pass is needed.
pub fn min_hops(
seeds: &Vector<bool>,
a: &Matrix<bool>,
max_hops: Option<u32>,
) -> Result<Vector<u64>> {
let n = square(a)?;
same_size(seeds, n, "seeds")?;
let idx = seeds.indices()?;
let dist = Vector::<u64>::from_tuples(n, &idx, &vec![0; idx.len()])?;
let frontier = dist.dup()?;
let mut hops = 0u32;
while frontier.nvals()? > 0 && max_hops.is_none_or(|h| hops < h) {
// frontier⟨¬dist, struct, replace⟩ = frontier MIN.PLUS A (bool edge = 1 hop)
grb!(GrB_vxm(
frontier.raw,
dist.raw,
null_mut(),
sys::GrB_MIN_PLUS_SEMIRING_UINT64,
frontier.raw,
a.raw,
sys::GrB_DESC_RSC
))?;
// dist⟨frontier, struct⟩ = frontier
grb!(GrB_Vector_assign(
dist.raw,
frontier.raw,
null_mut(),
frontier.raw,
sys::GrB_ALL,
n,
sys::GrB_DESC_S
))?;
hops += 1;
}
Ok(dist)
}
/// [`min_hops`] that STOPS after the first hop layer holding a node of
/// `targets` (or the seeds, when one is a target): every distance up to that
/// layer is the one [`min_hops`] computes, nothing beyond it is — what a
/// nearest-target path needs (the fewest hops to the set, then the hop
/// layers back), without walking everything the seeds reach.
pub fn min_hops_to(
seeds: &Vector<bool>,
a: &Matrix<bool>,
max_hops: Option<u32>,
targets: &[u64],
) -> Result<Vector<u64>> {
let n = square(a)?;
same_size(seeds, n, "seeds")?;
let idx = seeds.indices()?;
let dist = Vector::<u64>::from_tuples(n, &idx, &vec![0; idx.len()])?;
let hit = |v: &Vector<u64>| -> Result<bool> {
for t in targets {
if v.get(*t)?.is_some() {
return Ok(true);
}
}
Ok(false)
};
if hit(&dist)? {
return Ok(dist);
}
let frontier = dist.dup()?;
let mut hops = 0u32;
while frontier.nvals()? > 0 && max_hops.is_none_or(|h| hops < h) {
grb!(GrB_vxm(
frontier.raw,
dist.raw,
null_mut(),
sys::GrB_MIN_PLUS_SEMIRING_UINT64,
frontier.raw,
a.raw,
sys::GrB_DESC_RSC
))?;
grb!(GrB_Vector_assign(
dist.raw,
frontier.raw,
null_mut(),
frontier.raw,
sys::GrB_ALL,
n,
sys::GrB_DESC_S
))?;
hops += 1;
if hit(&frontier)? {
break;
}
}
Ok(dist)
}
/// For every node, the maximum of `vals` over all nodes it reaches (itself
/// included): "worst severity on any path". Fixpoint of
/// `s = max(s, A MAX.SECOND s)`; nodes reaching no valued node stay absent.
pub fn max_on_path(a: &Matrix<bool>, vals: &Vector<u64>) -> Result<Vector<u64>> {
let n = square(a)?;
same_size(vals, n, "vals")?;
let s = vals.dup()?;
let t = Vector::<u64>::new(n)?;
let improved = Vector::<bool>::new(n)?;
loop {
let before = s.nvals()?;
// t(i) = max over successors j of s(j)
grb!(GrB_mxv(
t.raw,
null_mut(),
null_mut(),
sys::GrB_MAX_SECOND_SEMIRING_UINT64,
a.raw,
s.raw,
sys::GrB_DESC_R
))?;
// improved⟨t⟩ = t > s ; new entries (only in t) count as improved
// even when t = 0, so mark them separately below.
grb!(GrB_Vector_eWiseMult_BinaryOp(
improved.raw,
null_mut(),
null_mut(),
sys::GrB_GT_UINT64,
t.raw,
s.raw,
sys::GrB_DESC_R
))?;
let mut any = false;
grb!(GrB_Vector_reduce_BOOL(
&mut any,
null_mut(),
sys::GrB_LOR_MONOID_BOOL,
improved.raw,
null_mut()
))?;
// s = max(s, t) (union: also adds nodes newly reaching a value)
grb!(GrB_Vector_eWiseAdd_BinaryOp(
s.raw,
null_mut(),
null_mut(),
sys::GrB_MAX_UINT64,
s.raw,
t.raw,
null_mut()
))?;
if !any && s.nvals()? == before {
break;
}
}
Ok(s)
}
/// The bottleneck walk — the (max, min) semiring over NODE weights (widest
/// path; "a chain is as strong as its weakest hop"). For every node a seed
/// reaches: the best, over all seed → node paths, of the SMALLEST weight on
/// the path (seed and node included). A node without a weight cannot be
/// passed and is absent from the result, as is an unreachable node.
///
/// Frontier fixpoint: `c = (f MAX.FIRST A) min w` proposes a value for every
/// successor of the frontier; a node stays in the frontier only where the
/// proposal beats what it had (or it had nothing). Values only ever rise and
/// are drawn from `weights`, so it terminates.
pub fn bottleneck(
seeds: &Vector<bool>,
a: &Matrix<bool>,
weights: &Vector<u64>,
) -> Result<Vector<u64>> {
let n = square(a)?;
same_size(seeds, n, "seeds")?;
same_size(weights, n, "weights")?;
// best⟨seeds⟩ = weights: a seed's own weight caps every path from it
let best = mask_vector(weights, seeds)?;
let frontier = best.dup()?;
let cand = Vector::<u64>::new(n)?;
let gt = Vector::<bool>::new(n)?;
while frontier.nvals()? > 0 {
// cand(j) = max over in-edges i → j of frontier(i)
grb!(GrB_vxm(
cand.raw,
null_mut(),
null_mut(),
sys::GxB_MAX_FIRST_UINT64,
frontier.raw,
a.raw,
sys::GrB_DESC_R
))?;
// cand = min(cand, weights): the hop is capped by the node entered
// (intersection: an unweighted node drops out)
grb!(GrB_Vector_eWiseMult_BinaryOp(
cand.raw,
null_mut(),
null_mut(),
sys::GrB_MIN_UINT64,
cand.raw,
weights.raw,
sys::GrB_DESC_R
))?;
// gt = cand > best, where both are present
grb!(GrB_Vector_eWiseMult_BinaryOp(
gt.raw,
null_mut(),
null_mut(),
sys::GrB_GT_UINT64,
cand.raw,
best.raw,
sys::GrB_DESC_R
))?;
// frontier⟨gt, replace⟩ = cand: the nodes whose value rose …
grb!(GrB_Vector_apply(
frontier.raw,
gt.raw,
null_mut(),
sys::GrB_IDENTITY_UINT64,
cand.raw,
sys::GrB_DESC_R
))?;
// … frontier⟨¬best, struct⟩ = cand: and the nodes reached for the
// first time (no replace: the entries above lie inside `best`)
grb!(GrB_Vector_apply(
frontier.raw,
best.raw,
null_mut(),
sys::GrB_IDENTITY_UINT64,
cand.raw,
sys::GrB_DESC_SC
))?;
// best = max(best, cand)
grb!(GrB_Vector_eWiseAdd_BinaryOp(
best.raw,
null_mut(),
null_mut(),
sys::GrB_MAX_UINT64,
best.raw,
cand.raw,
null_mut()
))?;
}
Ok(best)
}
/// `A` restricted to edges whose both ends are in `set` (dimensions kept):
/// `diag(set) ANY.SECOND A ANY.FIRST diag(set)`.
pub fn mask_by(a: &Matrix<bool>, set: &Vector<bool>) -> Result<Matrix<bool>> {
let n = square(a)?;
same_size(set, n, "set")?;
let mut d: sys::GrB_Matrix = null_mut();
grb!(GrB_Matrix_diag(&mut d, set.raw, 0))?;
let d = Matrix::<bool>::from_raw(d);
let rows = Matrix::<bool>::new(n, n)?;
let out = Matrix::<bool>::new(n, n)?;
grb!(GrB_mxm(
rows.raw,
null_mut(),
null_mut(),
sys::GxB_ANY_SECOND_BOOL,
d.raw,
a.raw,
null_mut()
))?;
grb!(GrB_mxm(
out.raw,
null_mut(),
null_mut(),
sys::GxB_ANY_FIRST_BOOL,
rows.raw,
d.raw,
null_mut()
))?;
Ok(out)
}
/// `v` restricted to the indices in `set` (`w⟨set, struct⟩ = v`).
pub fn mask_vector<T: crate::Scalar>(v: &Vector<T>, set: &Vector<bool>) -> Result<Vector<T>> {
let n = v.size()?;
same_size(set, n, "set")?;
let w = Vector::<T>::new(n)?;
let id = identity::<T>();
grb!(GrB_Vector_apply(
w.raw,
set.raw,
null_mut(),
id,
v.raw,
sys::GrB_DESC_S
))?;
Ok(w)
}
fn identity<T: crate::Scalar>() -> sys::GrB_UnaryOp {
let t = T::grb_type();
unsafe {
if t == sys::GrB_BOOL {
sys::GrB_IDENTITY_BOOL
} else if t == sys::GrB_UINT64 {
sys::GrB_IDENTITY_UINT64
} else {
sys::GrB_IDENTITY_FP64
}
}
}
/// Structural similarity `C = A·Aᵀ`: `C[i][j]` = number of successors `i` and
/// `j` share (diagonal = out-degree). PLUS.PAIR, which equals PLUS.TIMES on a
/// boolean pattern and skips reading values.
pub fn similarity(a: &Matrix<bool>) -> Result<Matrix<u64>> {
let n = a.nrows()?;
let c = Matrix::<u64>::new(n, n)?;
grb!(GrB_mxm(
c.raw,
null_mut(),
null_mut(),
sys::GxB_PLUS_PAIR_UINT64,
a.raw,
a.raw,
sys::GrB_DESC_T1
))?;
Ok(c)
}
/// One propagation step over a WEIGHTED adjacency: `y = A·x`, PLUS.TIMES on
/// FP64 — `y[i] = Σ_j A[i][j]·x[j]`. With `A` row-normalised it is the mean
/// of the neighbours' values; applied per column of a dense block it is the
/// propagation of FastRP (Chen et al., *Fast and Accurate Network Embeddings
/// via Very Sparse Random Projection*, CIKM 2019). A row with no neighbour has
/// no entry in `y` (read it as 0).
pub fn propagate(a: &Matrix<f64>, x: &Vector<f64>) -> Result<Vector<f64>> {
let n = a.nrows()?;
same_size(x, a.ncols()?, "x")?;
let y = Vector::<f64>::new(n)?;
grb!(GrB_mxv(
y.raw,
null_mut(),
null_mut(),
sys::GrB_PLUS_TIMES_SEMIRING_FP64,
a.raw,
x.raw,
null_mut()
))?;
Ok(y)
}
/// Weakly connected components: label = smallest node id in the component.
/// Label propagation over `A ∨ Aᵀ` with MIN.SECOND, plus one pointer jump
/// (`l = l[l]`) per round, until no label changes.
pub fn components(a: &Matrix<bool>) -> Result<Vector<u64>> {
let n = square(a)?;
let sym = Matrix::<bool>::new(n, n)?;
grb!(GrB_Matrix_eWiseAdd_BinaryOp(
sym.raw,
null_mut(),
null_mut(),
sys::GrB_LOR,
a.raw,
a.raw,
sys::GrB_DESC_T1
))?;
let ids: Vec<u64> = (0..n).collect();
let labels = Vector::<u64>::from_tuples(n, &ids, &ids)?;
let t = Vector::<u64>::new(n)?;
let changed = Vector::<bool>::new(n)?;
loop {
let mut any = false;
// t(i) = min label among neighbours
grb!(GrB_mxv(
t.raw,
null_mut(),
null_mut(),
sys::GxB_MIN_SECOND_UINT64,
sym.raw,
labels.raw,
sys::GrB_DESC_R
))?;
grb!(GrB_Vector_eWiseMult_BinaryOp(
changed.raw,
null_mut(),
null_mut(),
sys::GrB_LT_UINT64,
t.raw,
labels.raw,
sys::GrB_DESC_R
))?;
grb!(GrB_Vector_reduce_BOOL(
&mut any,
null_mut(),
sys::GrB_LOR_MONOID_BOOL,
changed.raw,
null_mut()
))?;
grb!(GrB_Vector_eWiseAdd_BinaryOp(
labels.raw,
null_mut(),
null_mut(),
sys::GrB_MIN_UINT64,
labels.raw,
t.raw,
null_mut()
))?;
if !any {
break;
}
// Pointer jump: labels(i) = labels(labels(i)). Labels only ever point
// to a smaller id in the same component, so this is safe and shortens
// the number of rounds from O(diameter) towards O(log n).
let (_, l) = labels.to_tuples()?;
let jumped = Vector::<u64>::new(n)?;
grb!(GrB_Vector_extract(
jumped.raw,
null_mut(),
null_mut(),
labels.raw,
l.as_ptr(),
n,
null_mut()
))?;
grb!(GrB_Vector_eWiseAdd_BinaryOp(
labels.raw,
null_mut(),
null_mut(),
sys::GrB_MIN_UINT64,
labels.raw,
jumped.raw,
null_mut()
))?;
}
Ok(labels)
}
/// Edges added and removed between two snapshots of the same node space:
/// `(added = Ay⟨¬Ax⟩, removed = Ax⟨¬Ay⟩)`.
pub fn diff(ax: &Matrix<bool>, ay: &Matrix<bool>) -> Result<(Matrix<bool>, Matrix<bool>)> {
let (n, m) = (ax.nrows()?, ax.ncols()?);
if (ay.nrows()?, ay.ncols()?) != (n, m) {
return Err(Error::invalid("diff: snapshots have different dimensions"));
}
let added = Matrix::<bool>::new(n, m)?;
let removed = Matrix::<bool>::new(n, m)?;
grb!(GrB_Matrix_apply(
added.raw,
ax.raw,
null_mut(),
sys::GrB_IDENTITY_BOOL,
ay.raw,
sys::GrB_DESC_RSC
))?;
grb!(GrB_Matrix_apply(
removed.raw,
ay.raw,
null_mut(),
sys::GrB_IDENTITY_BOOL,
ax.raw,
sys::GrB_DESC_RSC
))?;
Ok((added, removed))
}
/// `A ∨ B`: every edge of either graph (eWiseAdd with LOR). Used to walk
/// several edge kinds as one adjacency (e.g. CALLS ∨ USES_TYPE).
pub fn union(a: &Matrix<bool>, b: &Matrix<bool>) -> Result<Matrix<bool>> {
let (n, m) = (a.nrows()?, a.ncols()?);
if (b.nrows()?, b.ncols()?) != (n, m) {
return Err(Error::invalid("union: graphs have different dimensions"));
}
let c = Matrix::<bool>::new(n, m)?;
grb!(GrB_Matrix_eWiseAdd_BinaryOp(
c.raw,
null_mut(),
null_mut(),
sys::GrB_LOR,
a.raw,
b.raw,
null_mut()
))?;
Ok(c)
}
/// `A` restricted to the ROWS in `set` (out-edges of the set's nodes, to
/// anywhere): `diag(set) ANY.SECOND A`.
pub fn mask_rows(a: &Matrix<bool>, set: &Vector<bool>) -> Result<Matrix<bool>> {
let n = square(a)?;
same_size(set, n, "set")?;
let mut d: sys::GrB_Matrix = null_mut();
grb!(GrB_Matrix_diag(&mut d, set.raw, 0))?;
let d = Matrix::<bool>::from_raw(d);
let out = Matrix::<bool>::new(n, n)?;
grb!(GrB_mxm(
out.raw,
null_mut(),
null_mut(),
sys::GxB_ANY_SECOND_BOOL,
d.raw,
a.raw,
null_mut()
))?;
Ok(out)
}
/// Many independent single-seed walks at once: entry `i` of the result is
/// what [`reach`] from `seeds[i]` alone returns (seed included, ascending).
///
/// `batch` seeds at a time form the rows of ONE frontier matrix `F` (B × n,
/// row i = the frontier of walk i) and its visited matrix `V`; a step is
/// `F⟨¬V, struct, replace⟩ = F ANY.PAIR A` (one mxm instead of B vxm) and
/// `V⟨F, struct⟩ = true`. All rows start together, so hop `k` of the matrix
/// is hop `k` of every walk: `max_hops` bounds each row exactly as it bounds
/// a single walk. For a reverse walk pass the cached `Aᵀ`.
pub fn reach_each(seeds: &[u64], a: &Matrix<bool>, max_hops: Option<u32>, batch: usize) -> Result<Vec<Vec<u64>>> {
let sets: Vec<&[u64]> = seeds.iter().map(std::slice::from_ref).collect();
reach_sets(&sets, a, max_hops, batch)
}
/// [`reach_each`] with a seed SET per walk: entry `i` is what [`reach`]
/// from all of `sets[i]` returns (seeds included, ascending).
pub fn reach_sets(sets: &[&[u64]], a: &Matrix<bool>, max_hops: Option<u32>, batch: usize) -> Result<Vec<Vec<u64>>> {
let n = square(a)?;
let batch = batch.max(1);
let mut out: Vec<Vec<u64>> = Vec::with_capacity(sets.len());
for chunk in sets.chunks(batch) {
let b = chunk.len() as u64;
let rows: Vec<u64> = chunk.iter().enumerate().flat_map(|(i, s)| std::iter::repeat_n(i as u64, s.len())).collect();
let cols: Vec<u64> = chunk.iter().flat_map(|s| s.iter().copied()).collect();
let frontier = Matrix::<bool>::from_edges(b, n, &rows, &cols)?;
let visited = frontier.dup()?;
let mut hops = 0u32;
loop {
if max_hops.is_some_and(|h| hops >= h) || frontier.nvals()? == 0 {
break;
}
// frontier⟨¬visited, struct, replace⟩ = frontier ANY.PAIR A
grb!(GrB_mxm(
frontier.raw,
visited.raw,
null_mut(),
sys::GxB_ANY_PAIR_BOOL,
frontier.raw,
a.raw,
sys::GrB_DESC_RSC
))?;
// visited⟨frontier, struct⟩ = true
grb!(GrB_Matrix_assign_BOOL(
visited.raw,
frontier.raw,
null_mut(),
true,
sys::GrB_ALL,
b,
sys::GrB_ALL,
n,
sys::GrB_DESC_S
))?;
hops += 1;
}
let (i, j, _) = visited.to_tuples()?;
let mut per: Vec<Vec<u64>> = vec![Vec::new(); chunk.len()];
for (r, c) in i.into_iter().zip(j) {
per[r as usize].push(c);
}
for mut p in per {
p.sort_unstable();
out.push(p);
}
}
Ok(out)
}