radiate_rust/engines/schema/
subset.rs

1use rand::{rngs::ThreadRng, Rng};
2
3pub fn individual_indexes(
4    random: &mut ThreadRng,
5    index: usize,
6    size: usize,
7    order: usize,
8) -> Vec<usize> {
9    let mut sub_set = subset(size as usize, order as usize, random);
10    let mut i = 0;
11    while sub_set[i] < index as i32 && i < sub_set.len() - 1 {
12        i += 1;
13    }
14
15    sub_set[i] = index as i32;
16    sub_set.iter().map(|x| *x as usize).collect()
17}
18
19pub fn subset(n: usize, k: usize, random: &mut ThreadRng) -> Vec<i32> {
20    if k <= 0 {
21        panic!("Subset size smaller or equal zero: {}", k);
22    }
23
24    if n < k {
25        panic!("n smaller than k: {} < {}.", n, k);
26    }
27
28    let mut sub = vec![0; k as usize];
29    next(n as i32, &mut sub, random);
30    sub
31}
32
33fn next(num: i32, a: &mut Vec<i32>, random: &mut ThreadRng) {
34    let k = a.len() as i32;
35
36    if k == num {
37        for i in 0..k {
38            a[i as usize] = i;
39        }
40
41        return;
42    }
43
44    if k > num - k {
45        build_subset(num, a, random);
46        invert(num, a);
47    } else {
48        build_subset(num, a, random);
49    }
50}
51
52fn build_subset(n: i32, sub: &mut Vec<i32>, random: &mut ThreadRng) {
53    let k = sub.len() as i32;
54    check_subset(n, k);
55
56    if sub.len() as i32 == n {
57        for i in 0..k {
58            sub[i as usize] = i;
59        }
60        return;
61    }
62
63    for i in 0..k {
64        sub[i as usize] = i * n / k;
65    }
66
67    let mut l;
68    let mut ix;
69    for _ in 0..k {
70        loop {
71            ix = random.gen_range(1..n);
72            l = (ix * k - 1) / n;
73            if sub[l as usize] < ix {
74                break;
75            }
76        }
77
78        sub[l as usize] = sub[l as usize] + 1;
79    }
80
81    let mut m = 0;
82    let mut ip = 0;
83    let mut is_ = k;
84    for i in 0..k {
85        m = sub[i as usize];
86        sub[i as usize] = 0;
87
88        if m != i * n / k {
89            ip = ip + 1;
90            sub[ip as usize - 1] = m;
91        }
92    }
93
94    let ihi = ip;
95    for i in 1..=ihi {
96        ip = ihi + 1 - i;
97        l = 1 + (sub[ip as usize - 1] * k - 1) / n;
98        let ids = sub[ip as usize - 1] - (l - 1) * n / k;
99        sub[ip as usize - 1] = 0;
100        sub[is_ as usize - 1] = l;
101        is_ = is_ - ids;
102    }
103
104    let mut ir = 0;
105    let mut m0 = 0;
106    for ll in 1..=k {
107        l = k + 1 - ll;
108
109        if sub[l as usize - 1] != 0 {
110            ir = l;
111            m0 = 1 + (sub[l as usize - 1] - 1) * n / k;
112            m = sub[l as usize - 1] * n / k - m0 + 1;
113        }
114
115        ix = random.gen_range(m0..m0 + m - 1);
116
117        let mut i = l + 1;
118        while i <= ir && ix >= sub[i as usize - 1] {
119            ix = ix + 1;
120            sub[i as usize - 2] = sub[i as usize - 1];
121            i = i + 1;
122        }
123
124        sub[i as usize - 2] = ix;
125        m -= 1;
126    }
127}
128
129fn invert(n: i32, a: &mut Vec<i32>) {
130    let k = a.len() as i32;
131
132    let mut v = n - 1;
133    let mut j = n - k - 1;
134    let mut vi;
135
136    let mut ac = vec![0; k as usize];
137    ac.copy_from_slice(&a);
138
139    for i in (0..k).rev() {
140        loop {
141            vi = index_of(&ac, j, v);
142            if vi != -1 {
143                break;
144            }
145            v -= 1;
146            j = vi;
147        }
148
149        a[i as usize] = v;
150        v -= 1;
151    }
152}
153
154fn index_of(a: &Vec<i32>, start: i32, value: i32) -> i32 {
155    for i in (0..=start).rev() {
156        if a[i as usize] < value {
157            return -1;
158        } else if a[i as usize] == value {
159            return i;
160        }
161    }
162
163    -1
164}
165
166fn check_subset(n: i32, k: i32) {
167    if k <= 0 {
168        panic!("Subset size smaller or equal zero: {}", k);
169    }
170    if n < k {
171        panic!("n smaller than k: {} < {}.", n, k);
172    }
173}