pub mod arrayadapter;
pub mod safeadd;
pub use crate::arrayadapter::ArrayAdapter;
pub use crate::safeadd::SafeAdd;
use num_traits::{NumAssignOps, Signed, Zero};
#[derive(Debug, Copy, Clone)]
struct DistancePair<N> {
i: u32,
d: N,
}
#[derive(Debug)]
struct Rec<N> {
near: DistancePair<N>,
seco: DistancePair<N>,
}
#[inline]
fn initial_assignment<M, N>(mat: &M, med: &Vec<usize>, data: &mut Vec<Rec<N>>) -> N
where
N: NumAssignOps + Zero + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let n = mat.len();
let k = med.len();
let firstcenter = med[0];
let unassigned = k as u32;
let mut loss: N = N::zero();
for i in 0..n {
let mut cur = Rec::<N> {
near: DistancePair {
i: 0,
d: mat.get(i, firstcenter),
},
seco: DistancePair {
i: unassigned,
d: N::zero(),
},
};
for m in 1..k {
let d = mat.get(i, med[m]);
if d < cur.near.d || i == med[m] {
cur.seco = cur.near;
cur.near = DistancePair { i: m as u32, d: d };
} else if cur.seco.i == unassigned || d < cur.seco.d {
cur.seco = DistancePair { i: m as u32, d: d };
}
}
loss.safe_inc(cur.near.d);
data.push(cur);
}
return loss;
}
fn update_removal_loss<N>(data: &Vec<Rec<N>>, loss: &mut Vec<N>)
where
N: NumAssignOps + Signed + Copy + Zero + SafeAdd,
{
let n = data.len();
for i in 0..loss.len() {
loss[i] = N::zero();
}
for i in 0..n {
let rec = &data[i];
loss[rec.near.i as usize].safe_inc(rec.seco.d - rec.near.d);
}
}
#[inline]
fn find_min<N>(a: &Vec<N>) -> (usize, N)
where
N: PartialOrd + Copy + Zero,
{
let mut rk: usize = a.len();
let mut rv: N = N::zero();
for (ik, iv) in a.iter().enumerate() {
if ik == 0 || *iv < rv {
rk = ik;
rv = *iv;
}
}
return (rk, rv);
}
#[inline]
fn update_second_nearest<M, N>(
mat: &M,
med: &Vec<usize>,
n: usize,
b: usize,
o: usize,
djo: N,
) -> DistancePair<N>
where
N: NumAssignOps + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let mut s = DistancePair {
i: b as u32,
d: djo,
};
for i in 0..med.len() {
if i == n || i == b {
continue;
}
let dm = mat.get(o, med[i]);
if dm < s.d {
s = DistancePair { i: i as u32, d: dm };
}
}
return s;
}
#[inline]
fn find_best_swap<M, N>(mat: &M, removal_loss: &Vec<N>, data: &Vec<Rec<N>>, j: usize) -> (usize, N)
where
N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let n = mat.len();
let mut ploss = removal_loss.clone();
let mut acc = N::zero();
for o in 0..n {
let reco = &data[o];
let djo = mat.get(j, o);
if djo < reco.near.d {
acc.safe_inc(djo - reco.near.d);
ploss[reco.near.i as usize].safe_inc(reco.near.d - reco.seco.d);
} else if djo < reco.seco.d {
ploss[reco.near.i as usize].safe_inc(djo - reco.seco.d);
}
}
let (b, bloss) = find_min(&ploss);
return (b, bloss + acc);
}
#[cfg(feature = "assertions")]
fn debug_validate_assignment<M, N>(mat: &M, med: &Vec<usize>, data: &Vec<Rec<N>>)
where
N: NumAssignOps + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let n = mat.len();
for o in 0..n {
debug_assert!(
mat.get(o, med[data[o].near.i as usize]) == data[o].near.d,
"primary assignment inconsistent"
);
debug_assert!(
mat.get(o, med[data[o].seco.i as usize]) == data[o].seco.d,
"secondary assignment inconsistent"
);
debug_assert!(
data[o].near.d <= data[o].seco.d,
"nearest is farther than second nearest"
);
}
}
#[inline]
fn do_swap<M, N>(mat: &M, med: &mut Vec<usize>, data: &mut Vec<Rec<N>>, b: usize, j: usize) -> N
where
N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let n = mat.len();
med[b] = j;
let mut newloss = N::zero();
for o in 0..n {
let mut reco = &mut data[o];
if o == j {
if reco.near.i != b as u32 {
reco.seco = reco.near;
}
reco.near = DistancePair {
i: b as u32,
d: N::zero(),
};
continue;
}
let djo = mat.get(j, o);
if reco.near.i == b as u32 {
if djo < reco.seco.d {
reco.near = DistancePair {
i: b as u32,
d: djo,
};
} else {
reco.near = reco.seco;
reco.seco = update_second_nearest(mat, &med, reco.near.i as usize, b, o, djo);
}
} else {
if djo < reco.near.d {
reco.seco = reco.near;
reco.near = DistancePair {
i: b as u32,
d: djo,
};
} else if reco.seco.i == b as u32 {
reco.seco = update_second_nearest(mat, &med, reco.near.i as usize, b, o, djo);
} else if djo < reco.seco.d {
reco.seco = DistancePair {
i: b as u32,
d: djo,
};
}
}
newloss.safe_inc(reco.near.d);
}
#[cfg(feature = "assertions")]
debug_validate_assignment(&mat, &med, &data);
return newloss;
}
#[cfg(feature = "rand")]
pub fn random_initialization(n: usize, k: usize, mut rng: &mut impl rand::Rng) -> Vec<usize> {
return rand::seq::index::sample(&mut rng, n, k).into_vec();
}
pub fn fasterpam<M, N>(
mat: &M,
mut med: &mut Vec<usize>,
maxiter: usize,
) -> (N, usize, usize, Vec<usize>)
where
N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd,
M: ArrayAdapter<N>,
{
let n = mat.len();
let k = med.len();
assert!(mat.is_square(), "Dissimilarity matrix is not square");
assert!(n <= u32::MAX as usize, "N is too large");
assert!(k > 0 && k < u32::MAX as usize, "invalid N");
assert!(k <= n, "k must be at most N");
let mut data = Vec::<Rec<N>>::with_capacity(n);
let mut loss = initial_assignment(mat, &med, &mut data);
#[cfg(feature = "assertions")]
debug_validate_assignment(&mat, &med, &data);
let mut removal_loss = vec![N::zero(); k];
let mut lastswap = n;
let mut numswaps = 0;
let mut iter = 0;
while iter < maxiter {
iter += 1;
let swaps_before = numswaps;
update_removal_loss(&data, &mut removal_loss);
for j in 0..n {
if j == lastswap {
break;
}
if j == data[j].near.i as usize {
continue;
}
let (b, change) = find_best_swap(mat, &removal_loss, &data, j);
if change >= N::zero() {
continue;
}
numswaps += 1;
lastswap = j;
let newloss = do_swap(mat, &mut med, &mut data, b, j);
if newloss >= loss {
break;
}
loss = newloss;
update_removal_loss(&data, &mut removal_loss);
}
if numswaps == swaps_before {
break;
}
}
let assignment = data.iter().map(|x| x.near.i as usize).collect();
return (loss, numswaps, iter, assignment);
}
#[cfg(test)]
mod tests {
use crate::arrayadapter::LowerTriangle;
use crate::fasterpam;
#[test]
fn basic() {
let data = LowerTriangle {
n: 5,
data: vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 1],
};
let mut meds = vec![0, 1];
let (loss, numswap, numiter, assignment) = fasterpam(&data, &mut meds, 10);
assert_eq!(loss, 4, "loss not as expected");
assert_eq!(numswap, 2, "swaps not as expected");
assert_eq!(numiter, 2, "iterations not as expected");
let expected = vec![0, 0, 0, 1, 1];
assert!(
assignment.iter().zip(expected.iter()).all(|(a, b)| a == b),
"cluster assignment not as expected"
);
}
}