pub mod safeadd;
pub use crate::safeadd::SafeAdd;
use ndarray::Array2;
use num_traits::{NumAssignOps, Signed, Zero};
#[derive(Debug, Copy, Clone)]
struct DistancePair<N> {
i: u32,
d: N,
}
#[derive(Debug)]
struct PointInformation<N> {
near: DistancePair<N>,
seco: DistancePair<N>,
}
fn initial_assignment<N: NumAssignOps + Zero + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
med: &Vec<usize>,
data: &mut Vec<PointInformation<N>>,
) -> N {
let n = mat.shape()[0];
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 = PointInformation::<N> {
near: DistancePair {
i: 0,
d: mat[[i, firstcenter]],
},
seco: DistancePair {
i: unassigned,
d: N::zero(),
},
};
for m in 1..k {
let d = mat[[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: NumAssignOps + Signed + Copy + Zero + SafeAdd>(
data: &Vec<PointInformation<N>>,
loss: &mut Vec<N>,
) {
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: PartialOrd + Copy + Zero>(a: &Vec<N>) -> (usize, N) {
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);
}
fn update_second_nearest<N: NumAssignOps + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
med: &Vec<usize>,
n: usize,
b: usize,
o: usize,
djo: N,
) -> DistancePair<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[[o, med[i]]];
if dm < s.d {
s = DistancePair { i: i as u32, d: dm };
}
}
return s;
}
fn find_best_swap<N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
removal_loss: &Vec<N>,
data: &Vec<PointInformation<N>>,
j: usize,
) -> (usize, N) {
let n = mat.shape()[0];
let mut ploss = removal_loss.clone();
let mut acc = N::zero();
for o in 0..n {
let reco = &data[o];
let djo = mat[[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<N: NumAssignOps + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
med: &Vec<usize>,
data: &Vec<PointInformation<N>>,
) {
let n = mat.shape()[0];
for o in 0..n {
debug_assert!(
mat[[o, med[data[o].near.i as usize]]] == data[o].near.d,
"primary assignment inconsistent"
);
debug_assert!(
mat[[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);
}
}
fn do_swap<N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
med: &mut Vec<usize>,
data: &mut Vec<PointInformation<N>>,
b: usize,
j: usize,
) -> N {
let n = mat.shape()[0];
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[[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;
}
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<N: NumAssignOps + Signed + Zero + PartialOrd + Copy + SafeAdd>(
mat: &Array2<N>,
mut med: &mut Vec<usize>,
maxiter: usize,
) -> (N, usize, usize, Vec<usize>) {
let n = mat.shape()[0];
let k = med.len();
assert_eq!(mat.shape()[1], n);
assert!(k <= n);
let mut data = Vec::<PointInformation<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);
}