001/*
002 * Java Genetic Algorithm Library (jenetics-7.2.0).
003 * Copyright (c) 2007-2023 Franz Wilhelmstötter
004 *
005 * Licensed under the Apache License, Version 2.0 (the "License");
006 * you may not use this file except in compliance with the License.
007 * You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 *
017 * Author:
018 *    Franz Wilhelmstötter (franz.wilhelmstoetter@gmail.com)
019 */
020package io.jenetics.internal.math;
021
022import static java.lang.String.format;
023import static java.util.Objects.requireNonNull;
024import static io.jenetics.internal.math.Basics.isMultiplicationSave;
025
026import java.util.Arrays;
027import java.util.random.RandomGenerator;
028
029/**
030 * This class creates random subsets of size {@code k}  from a set of {@code n}
031 * elements. Implementation of the {@code RANKSB} algorithm described by
032 * <em>Albert Nijenhuis</em> and <em>Herbert Wilf</em> in <b>Combinatorial
033 * Algorithms for Computers and Calculators</b>
034 * <p>
035 *  Reference:<em>
036 *      Albert Nijenhuis, Herbert Wilf,
037 *      Combinatorial Algorithms for Computers and Calculators,
038 *      Second Edition,
039 *      Academic Press, 1978,
040 *      ISBN: 0-12-519260-6,
041 *      LC: QA164.N54.
042 *      Page: 42</em>
043 *
044 * @author <a href="mailto:franz.wilhelmstoetter@gmail.com">Franz Wilhelmstötter</a>
045 * @version 7.0
046 * @since 4.0
047 */
048public final class Subset {
049        private Subset() {}
050
051        /**
052         * Selects a random subset of size {@code k} from a set of size {@code n}.
053         *
054         * @param n the size of the set.
055         * @param k the size of the subset.
056         * @param rnd the random number generator used.
057         * @throws NullPointerException if {@code random} is {@code null}.
058         * @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
059         *         {@code n*k} will cause an integer overflow.
060         * @return the a-set array for the given parameter. The returned sub-set
061         *         array is sorted in increasing order.
062         */
063        public static int[] next(final RandomGenerator rnd, final int n, final int k) {
064                requireNonNull(rnd);
065
066                final var subset = new int[k];
067                next(rnd, n, subset);
068                return subset;
069        }
070
071        private static void next(
072                final RandomGenerator rnd,
073                final int n,
074                final int[] a
075        ) {
076                requireNonNull(rnd, "Random");
077                requireNonNull(a, "Sub set array");
078
079                final int k = a.length;
080                checkSubSet(n, k);
081
082                // Early return if k == n.
083                if (k == n) {
084                        for (int i = 0; i < k; ++i) {
085                                a[i] = i;
086                        }
087                        return;
088                }
089
090                // Calculate the 'inverse' subset if k > n - k.
091                if (k > n - k) {
092                        subset0(rnd, n, n - k, a);
093                        invert(n, a);
094                } else {
095                        subset0(rnd, n, k, a);
096                }
097        }
098
099        /*
100         * "Inverts" the given subset array `a`. The first n - k elements represent
101         * the set, which must not be part of the "inverted" subset. This is done by
102         * filling the array from the back, starting with the highest possible element,
103         * which is not part of the "forbidden" subset elements. The result is a
104         * subset array, filled with elements, which is not part of the original
105         * "forbidden" subset.
106         */
107        private static void invert(
108                final int n,
109                final int[] a
110        ) {
111                final int k = a.length;
112
113                int v = n - 1;
114                int j = n - k - 1;
115                int vi;
116
117                final int[] ac = Arrays.copyOfRange(a, 0, n - k);
118                for (int i = k; --i >= 0;) {
119                        while ((vi = indexOf(ac, j, v)) != -1) {
120                                --v;
121                                j = vi;
122                        }
123
124                        a[i] = v--;
125                }
126        }
127
128        private static int indexOf(final int[] a, final int start, final int value) {
129                for (int i = start; i >= 0; --i) {
130                        if (a[i] < value) {
131                                return -1;
132                        } else if (a[i] == value) {
133                                return i;
134                        }
135                }
136
137                return -1;
138        }
139
140        // The actual implementation of the `RANKSB` algorithm.
141        private static void subset0(
142                final RandomGenerator random,
143                final int n,
144                final int k,
145                final int[] a
146        ) {
147                assert k <= a.length;
148                assert k <= n - k;
149
150                // Early return if k = 1.
151                if (k == 1) {
152                        a[0] = random.nextInt(n);
153                        return;
154                }
155
156                // (A): Initialize a[i] to "zero" point for bin Ri.
157                for (int i = 0; i < k; ++i) {
158                        a[i] = (i*n)/k;
159                }
160
161                // (B)
162                int l, x;
163                for (int c = 0; c < k; ++c) {
164                        do {
165                                // Choose random x;
166                                x = 1 + random.nextInt(n);
167
168                                // determine range Rl;
169                                l = (x*k - 1)/n;
170                        } while (a[l] >= x); // accept or reject.
171
172                        ++a[l];
173                }
174                int s = k;
175
176                // (C) Move a[i] of nonempty bins to the left.
177                int m = 0, p = 0;
178                for (int i = 0; i < k; ++i) {
179                        if (a[i] == (i*n)/k) {
180                                a[i] = 0;
181                        } else {
182                                ++p;
183                                m = a[i];
184                                a[i] = 0;
185                                a[p - 1] = m;
186                        }
187                }
188
189                // (D) Determine l, set up space for Bl.
190                for (; p > 0; --p) {
191                        l = 1 + (a[p - 1]*k - 1)/n;
192                        final int ds = a[p - 1] - ((l - 1)*n)/k;
193                        a[p - 1] = 0;
194                        a[s - 1] = l;
195                        s -= ds;
196                }
197
198                // (E) If a[l] != 0, a new bin is to be processed.
199                int r = 0, m0 = 0;
200                for (int ll = 1; ll <= k; ++ll) {
201                        l = k + 1 - ll;
202
203                        if (a[l - 1] != 0) {
204                                r = l;
205                                m0 = 1 + ((a[l - 1] - 1)*n)/k;
206                                m = (a[l - 1]*n)/k - m0 + 1;
207                        }
208
209                        // (F) Choose a random x.
210                        x = m0 + random.nextInt(m);
211                        int i = l + 1;
212
213                        // (G) Check x against previously entered elements in bin;
214                        //     increment x as it jumps over elements <= x.
215                        while (i <= r && x >= a[i - 1]) {
216                                ++x;
217                                a[i - 2] = a[i - 1];
218                                ++i;
219                        }
220
221                        a[i - 2] = x;
222                        --m;
223                }
224
225                // Convert to zero-based indexed arrays.
226                for (int i = 0; i < k; ++i) {
227                        --a[i];
228                }
229        }
230
231
232        private static void checkSubSet(final int n, final int k) {
233                if (k <= 0) {
234                        throw new IllegalArgumentException(format(
235                                "Subset size smaller or equal zero: %s", k
236                        ));
237                }
238                if (n < k) {
239                        throw new IllegalArgumentException(format(
240                                "n smaller than k: %s < %s.", n, k
241                        ));
242                }
243                if (!isMultiplicationSave(n, k)) {
244                        throw new IllegalArgumentException(format(
245                                "n*sub.length > Integer.MAX_VALUE (%s*%s = %s > %s)",
246                                n, k, (long)n*(long)k, Integer.MAX_VALUE
247                        ));
248                }
249        }
250
251}