1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#![lower::apply(algebraic)]
use atools::prelude::*;
use lower::algebraic::math;
pub trait Closest<const N: usize> {
    fn closest(&self, color: [f32; N]) -> (f32, [f32; N], u32);
    fn best(&self, color: [f32; N]) -> [f32; N] {
        self.closest(color).1
    }
    fn nearest(&self, color: [f32; N]) -> u32 {
        self.closest(color).2
    }
    fn space(&self) -> f32;
}

#[inline(always)]
fn euclidean_distance<const N: usize>(f: [f32; N], with: [f32; N]) -> f32 {
    f.asub(with)
        .map(|x| x * x)
        .into_iter()
        .fold(0.0, |acc, x| acc + x)
}

#[inline(always)]
fn minwby<T: Copy, U: PartialOrd>(max: T, x: impl Iterator<Item = T>, extractor: fn(T) -> U) -> T {
    x.fold(max, |acc, x| {
        if extractor(acc) > extractor(x) {
            x
        } else {
            acc
        }
    })
}

impl<'a, const N: usize> Closest<N> for super::pal<'a, N> {
    /// o(nn)
    fn closest(&self, color: [f32; N]) -> (f32, [f32; N], u32) {
        minwby(
            (f32::MAX, [0.0; N], 0),
            (0..)
                .zip(&**self)
                .map(|(i, &x)| (euclidean_distance(x, color), x, i)),
            |x| x.0,
        )
    }

    /// o(nn)
    #[cold]
    fn space(&self) -> f32 {
        self.iter()
            .enumerate()
            .map(|(i, &x)| {
                minwby(
                    f32::MAX,
                    self.iter()
                        .enumerate()
                        .filter(|&(j, _)| j != i)
                        .map(|(_, &y)| euclidean_distance(y, x)),
                    std::convert::identity,
                )
                .sqrt()
            })
            .fold(0.0, |x, y| x + y)
            / self.len() as f32
    }
}