xicor 0.1.0

An implementation of Sourav Chatterjee's xi-correlation coefficient
Documentation 
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

use ordered_float::OrderedFloat;
use num_traits::float::FloatCore;



/// Calculate the normalised xi-correlation of two floating-point sequences.
///
/// See [`xicor_norm`] for details of normalisation.
///
/// # Example
///
/// ```
/// use xicor::xicorf_norm;
///
/// let x: Vec<f32> = (0..47).map(|i| i as f32).collect();
/// let y: Vec<f32> = x.iter().map(|x| x*x).collect();
/// let xi = xicorf_norm(&x, &y);
///
/// assert_eq!(xi, 1.);
/// ```
///
/// Note that this is exactly the same data used in the example for [`xicorf`],
/// but here the result is actually 1.
pub fn xicorf_norm<F: FloatCore>(x: &[F], y: &[F]) -> f64 {
    let n = x.len() as f64;
    let lim = (n-2.)/(n+1.);

    xicorf(x, y)/lim
}

/// Calculate the normalised xi-correlation of two sequences whose values are
/// orderable (they implement [`Ord`]).
///
/// One particular quirk of Chatterjee's Xi is that its maximum value is only 1
/// in the limit as the number of points goes to infinity. For a finite number
/// of points `n`, its upper limit is `(n-2)/(n+1)`, _even with perfect
/// correlation_. This can be substantially below 1 for small datasets.
/// Therefore, this function is provided which normalises by that limit, such
/// that it will always return 1 for perfectly correlated data.
///
/// # Example
///
/// ```
/// use xicor::xicor_norm;
///
/// let x: Vec<u32> = (0..47).collect();
/// let y: Vec<u32> = x.iter().map(|x| x*x).collect();
/// let xi = xicor_norm(&x, &y);
///
/// assert_eq!(xi, 1.);
/// ```
///
/// Note that this is exactly the same data used in the example for [`xicor`],
/// but here the result is actually 1.
pub fn xicor_norm<T: Ord + Copy>(x: &[T], y: &[T]) -> f64 {
    let n = x.len() as f64;
    let lim = (n-2.)/(n+1.);

    xicor(x, y)/lim
}

/// Calculate the xi-correlation of two floating-point sequences.
///
/// This is a thin wrapper around [`xicor`] that transmutes slices of floats
/// into slices of [`OrderedFloat`], which implements the necessary [`Ord`].
///
/// # Example
///
/// ```
/// use xicor::xicorf;
///
/// let x: Vec<f32> = (0..47).map(|i| i as f32).collect();
/// let y: Vec<f32> = x.iter().map(|x| x*x).collect();
/// let xi = xicorf(&x, &y);
///
/// assert_eq!(xi, 0.9375);
/// ```
pub fn xicorf<F: FloatCore>(x: &[F], y: &[F]) -> f64 {
    // This is safe because OrderedFloat has transparent representation
    let x: &[OrderedFloat<F>] = unsafe { std::mem::transmute(x) };
    let y: &[OrderedFloat<F>] = unsafe { std::mem::transmute(y) };

    xicor(x, y)
}

/// Calculate the xi-correlation of two sequences whose values are orderable
/// (they implement [`Ord`]).
///
/// # Example
///
/// ```
/// use xicor::xicor;
///
/// let x: Vec<u32> = (0..47).collect();
/// let y: Vec<u32> = x.iter().map(|x| x*x).collect();
/// let xi = xicor(&x, &y);
///
/// assert_eq!(xi, 0.9375);
/// ```
pub fn xicor<T: Ord + Copy>(x: &[T], y: &[T]) -> f64 {
    assert!(x.len() == y.len(), "x and y must have the same length");

    let idcs = argsort(x);
    let y_ord = permute(y, &idcs);

    let idcs = argsort(&y_ord);
    let y_ascending = permute(&y_ord, &idcs);
    let r_ascending = cumulative_lte(&y_ascending);
    let l_ascending = cumulative_gte(&y_ascending);
    let mut rs = vec![0.; x.len()];
    let mut ls = vec![0.; x.len()];

    for ((i, r), l) in idcs.into_iter().zip(r_ascending).zip(l_ascending) {
        rs[i] = r as f64;
        ls[i] = l as f64;
    }

    let rsum = rs.windows(2)
        .map(|win| (win[0]-win[1]).abs())
        .sum::<f64>();

    let n = x.len() as f64;
    let lsum = ls.into_iter()
        .map(|l| l*(n-l))
        .sum::<f64>();

    1.-n*rsum/(2.*lsum)
}

// Return the indices that would sort the given array. That is, if you map the
// returned sequence of indices i -> arr[i], the resulting sequence is sorted.
pub(super) fn argsort<T: Ord>(arr: &[T]) -> Vec<usize> {
    let mut idcs: Vec<usize> = (0..arr.len()).collect();

    idcs.sort_unstable_by_key(|&i| &arr[i]);
    idcs
}

// Permute the given array such that arr[n] ends up at idcs[n].
pub(super) fn permute<T: Copy>(arr: &[T], idcs: &[usize]) -> Vec<T> {
    idcs.iter()
        .map(|&i| arr[i])
        .collect()
}

// For every element in the array, count how many elements are less than or
// equal to it. The array should be sorted before it is passed in.
pub(super) fn cumulative_lte<T: PartialEq<T> + Copy>(arr: &[T]) -> Vec<usize> {
    let mut counts: Vec<usize> = (1..=arr.len()).collect();

    for i in (0..arr.len()-1).rev() {
        if arr[i] == arr[i+1] { counts[i] = counts[i+1]; }
    }

    counts
}

// For every element in the array, count how many elements are greater than or
// equal to it. The array should be sorted before it is passed in.
pub(super) fn cumulative_gte<T: PartialEq<T> + Copy>(arr: &[T]) -> Vec<usize> {
    let mut counts: Vec<usize> = (1..=arr.len()).rev().collect();

    for i in 0..arr.len()-1 {
        if arr[i+1] == arr[i] { counts[i+1] = counts[i]; }
    }

    counts
}