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
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

use num_traits::{NumOps, Zero};

pub(crate) fn det_correlation_matrix<T: Clone + NumOps + Zero>(mat: &[Vec<T>]) -> T {
    let dim = mat[0].len();
    let num = mat.len();
    let mut cor = Vec::new();
    for i in 0..num {
        let mut column = Vec::new();
        for j in 0..num {
            let mut c = T::zero();
            for k in 0..dim {
                c = c + mat[i][k].clone() * mat[j][k].clone();
            }
            column.push(c);
        }
        cor.push(column);
    }
    det(&cor)
}

pub(crate) fn dot<T: Clone + NumOps + Zero>(x: &[T], y: &[T]) -> T {
    x.iter()
        .zip(y.iter())
        .map(|(a, b)| a.clone() * b.clone())
        .fold(T::zero(), |sum, c| sum + c)
}

pub(crate) fn sub<T: Clone + NumOps + Zero>(x: &[T], y: &[T]) -> Vec<T> {
    x.iter()
        .zip(y.iter())
        .map(|(a, b)| a.clone() - b.clone())
        .collect()
}

pub(crate) fn min_max_index_each_axis<T: Clone + Zero + PartialOrd>(
    points: &[Vec<T>],
) -> Vec<(usize, usize)> {
    let dim = points[0].len();
    let mut min_index = vec![0; dim];
    let mut max_index = vec![0; dim];
    let mut min = vec![T::zero(); dim];
    let mut max = vec![T::zero(); dim];
    for (index, point) in points.iter().enumerate() {
        for (j, element) in point.iter().enumerate() {
            if index == 0 || *element < min[j] {
                min[j] = element.clone();
                min_index[j] = index;
            }
            if index == 0 || *element > max[j] {
                max[j] = element.clone();
                max_index[j] = index;
            }
        }
    }
    min_index.into_iter().zip(max_index.into_iter()).collect()
}

pub(crate) fn is_same_dimension<T>(points: &[Vec<T>]) -> bool {
    if points.len() == 0 {
        return true;
    }
    let dim = points[0].len();
    if points.iter().skip(1).find(|p| p.len() != dim).is_some() {
        false
    } else {
        true
    }
}

pub(crate) fn det<T: NumOps + Clone>(m: &[Vec<T>]) -> T {
    let row_dim = m.len();
    let column_dim = m[0].len();
    assert_eq!(row_dim, column_dim);
    match column_dim {
        1 => m[0][0].clone(),
        2 => det_2x2(m),
        3 => det_3x3(m),
        4 => det_4x4(m),
        _ => {
            unimplemented!("matrix size should be less than 4 dim to calcurate determinant for now")
        }
    }
}

fn det_2x2<T: NumOps + Clone>(m: &[Vec<T>]) -> T {
    m[0][0].clone() * m[1][1].clone() - m[0][1].clone() * m[1][0].clone()
}

#[rustfmt::skip]
fn det_3x3<T: NumOps+Clone>(m: &[Vec<T>]) -> T {
    m[0][0].clone() * (m[1][1].clone() * m[2][2].clone() - m[1][2].clone() * m[2][1].clone())
  - m[1][0].clone() * (m[0][1].clone() * m[2][2].clone() - m[0][2].clone() * m[2][1].clone())
  + m[2][0].clone() * (m[0][1].clone() * m[1][2].clone() - m[0][2].clone() * m[1][1].clone())
}

#[rustfmt::skip]
fn det_4x4<T: NumOps + Clone>(m: &[Vec<T>]) -> T {
    m[0][0].clone()
        * (   m[1][1].clone() * m[2][2].clone() * m[3][3].clone()
            + m[1][2].clone() * m[2][3].clone() * m[3][1].clone()
            + m[1][3].clone() * m[2][1].clone() * m[3][2].clone()
            - m[1][3].clone() * m[2][2].clone() * m[3][1].clone()
            - m[1][2].clone() * m[2][1].clone() * m[3][3].clone()
            - m[1][1].clone() * m[2][3].clone() * m[3][2].clone())
  - m[1][0].clone()
        * (   m[0][1].clone() * m[2][2].clone() * m[3][3].clone()
            + m[0][2].clone() * m[2][3].clone() * m[3][1].clone()
            + m[0][3].clone() * m[2][1].clone() * m[3][2].clone()
            - m[0][3].clone() * m[2][2].clone() * m[3][1].clone()
            - m[0][2].clone() * m[2][1].clone() * m[3][3].clone()
            - m[0][1].clone() * m[2][3].clone() * m[3][2].clone())
  + m[2][0].clone()
        * (   m[0][1].clone() * m[1][2].clone() * m[3][3].clone()
            + m[0][2].clone() * m[1][3].clone() * m[3][1].clone()
            + m[0][3].clone() * m[1][1].clone() * m[3][2].clone()
            - m[0][3].clone() * m[1][2].clone() * m[3][1].clone()
            - m[0][2].clone() * m[1][1].clone() * m[3][3].clone()
            - m[0][1].clone() * m[1][3].clone() * m[3][2].clone())
  - m[3][0].clone()
        * (   m[0][1].clone() * m[1][2].clone() * m[2][3].clone()
            + m[0][2].clone() * m[1][3].clone() * m[2][1].clone()
            + m[0][3].clone() * m[1][1].clone() * m[2][2].clone()
            - m[0][3].clone() * m[1][2].clone() * m[2][1].clone()
            - m[0][2].clone() * m[1][1].clone() * m[2][3].clone()
            - m[0][1].clone() * m[1][3].clone() * m[2][2].clone())
}

#[test]
fn det_4x4_test() {
    let r1 = vec![1., 2., 3., 4.];
    let r2 = vec![5., 6., 7., 8.];
    let r3 = vec![9., 10., 11., 12.];
    let r4 = vec![13., 14., 15., 16.];
    assert_eq!(det(&[r1, r2, r3, r4]), 0.);
    let r1 = vec![1., 2., 3., 1.];
    let r2 = vec![6., 6., 7., 1.];
    let r3 = vec![9., 11., 11., 1.];
    let r4 = vec![15., 14., 15., 1.];
    assert_eq!(det(&[r1, r2, r3, r4]), -4.);
    let r1 = vec![1., 1., 1., 1.];
    let r2 = vec![1., 2., 2., 1.];
    let r3 = vec![1., 2., 3., 1.];
    let r4 = vec![1., 1., 1., 2.];
    assert_eq!(det(&[r1, r2, r3, r4]), 1.);
}