imgal 0.2.0

A fast and open-source scientific image processing and algorithm library.
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

use std::cmp::Ordering;

use crate::error::ImgalError;
use crate::traits::numeric::AsNumeric;

/// Sort 1-dimensional arrays of values and their associated weights.
///
/// # Description
///
/// Performs a bottom up merge sort on the input 1-dimensional data array along
/// with it's associated weights. Both the `data` and `weights` arrays are
/// _mutated_ during the sorting. The output of this function is a weighted
/// inversion count.
///
/// # Arguments
///
/// * `data`: A 1-dimensional array/slice of numbers of the same length as
///   `weights`.
/// * `weights`: A 1-dimensional array/slice of weights of the same length as
///   `data`.
///
/// # Returns
///
/// * `OK(f64)`: The number of swaps needed to sort the input array.
/// * `Err(ImgalError)`: If the `data.len() != weights.len()`.
///
/// # Reference
///
/// <https://doi.org/10.1109/TIP.2019.2909194>
pub fn weighted_merge_sort_mut<T>(data: &mut [T], weights: &mut [f64]) -> Result<f64, ImgalError>
where
    T: AsNumeric,
{
    // validate the input arrays are the same length
    let dl = data.len();
    let wl = weights.len();
    if dl != wl {
        return Err(ImgalError::MismatchedArrayLengths {
            a_arr_name: "data",
            a_arr_len: dl,
            b_arr_name: "weights",
            b_arr_len: wl,
        });
    };

    // create sort parameters and working buffers, start weighted bottom-up
    // merge sort
    let mut swap = 0.0;
    let mut swap_temp: f64;
    let mut step: usize = 1;
    let mut left: usize;
    let mut right: usize;
    let mut end: usize;
    let mut k: usize;
    let mut data_buf = vec![T::default(); dl];
    let mut weights_buf = vec![0.0; dl];
    let mut cum_weights_buf = vec![0.0; dl];
    while step < dl {
        left = 0;
        k = 0;
        let mut cw_acc = weights[0];
        cum_weights_buf[0] = weights[0];
        cum_weights_buf
            .iter_mut()
            .zip(weights.iter())
            .skip(1)
            .for_each(|(cw, w)| {
                *cw = cw_acc + w;
                cw_acc = *cw;
            });

        loop {
            right = left + step;
            end = right + step;
            if end > dl {
                if right > dl {
                    break;
                }
                end = dl;
            }
            let mut l = left;
            let mut r = right;
            while l < right && r < end {
                match data[l].partial_cmp(&data[r]) {
                    Some(Ordering::Greater) => {
                        if l == 0 {
                            swap_temp = weights[r] * cum_weights_buf[right - 1];
                        } else {
                            swap_temp =
                                weights[r] * (cum_weights_buf[right - 1] - cum_weights_buf[l - 1]);
                        }
                        swap += swap_temp;
                        data_buf[k] = data[r];
                        weights_buf[k] = weights[r];
                        k += 1;
                        r += 1;
                    }
                    _ => {
                        data_buf[k] = data[l];
                        weights_buf[k] = weights[l];
                        k += 1;
                        l += 1;
                    }
                }
            }
            if l < right {
                while l < right {
                    data_buf[k] = data[l];
                    weights_buf[k] = weights[l];
                    k += 1;
                    l += 1;
                }
            } else {
                while r < end {
                    data_buf[k] = data[r];
                    weights_buf[k] = weights[r];
                    k += 1;
                    r += 1;
                }
            }
            left = end;
        }

        // copy any unmerged tail, if array size is not a power of 2
        if k < dl {
            while k < dl {
                data_buf[k] = data[k];
                weights_buf[k] = weights[k];
                k += 1;
            }
        }

        // prepare for the next step, copy merged results back source
        data.clone_from_slice(&data_buf);
        weights.clone_from_slice(&weights_buf);

        // double the run size, continue
        step *= 2;
    }

    Ok(swap)
}