001/*
002 * Java Genetic Algorithm Library (jenetics-7.1.0).
003 * Copyright (c) 2007-2022 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         * @see #next(int, int[], RandomGenerator)
055         *
056         * @param n the size of the set.
057         * @param k the size of the subset.
058         * @param rnd the random number generator used.
059         * @throws NullPointerException if {@code random} is {@code null}.
060         * @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
061         *         {@code n*k} will cause an integer overflow.
062         * @return the a-set array for the given parameter. The returned sub-set
063         *         array is sorted in increasing order.
064         */
065        public static int[] next(final int n, final int k, final RandomGenerator rnd) {
066                final var subset = new int[k];
067                next(n, subset, rnd);
068                return subset;
069        }
070
071        private static void next(
072                final int n,
073                final int[] a,
074                final RandomGenerator rnd
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(n, n - k, a, rnd);
093                        invert(n, k, a);
094                } else {
095                        subset0(n, k, a, rnd);
096                }
097        }
098
099        /*
100         * "Inverts" the given subset array `a`. The first n - k elements represents
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 where not part of the original
105         * "forbidden" subset.
106         */
107        private static void invert(
108                final int n,
109                final int k,
110                final int[] a
111        ) {
112                assert a.length == k;
113
114                int v = n - 1;
115                int j = n - k - 1;
116                int vi;
117
118                final int[] ac = Arrays.copyOfRange(a, 0, n - k);
119                for (int i = k; --i >= 0;) {
120                        while ((vi = indexOf(ac, j, v)) != -1) {
121                                --v;
122                                j = vi;
123                        }
124
125                        a[i] = v--;
126                }
127        }
128
129        private static int indexOf(final int[] a, final int start, final int value) {
130                for (int i = start; i >= 0; --i) {
131                        if (a[i] < value) {
132                                return -1;
133                        } else if (a[i] == value) {
134                                return i;
135                        }
136                }
137
138                return -1;
139        }
140
141        // The actual implementation of the `RANKSB` algorithm.
142        private static void subset0(
143                final int n,
144                final int k,
145                final int[] a,
146                final RandomGenerator random
147        ) {
148                assert k <= a.length;
149                assert k <= n - k;
150
151                // Early return if k = 1.
152                if (k == 1) {
153                        a[0] = random.nextInt(n);
154                        return;
155                }
156
157                // (A): Initialize a[i] to "zero" point for bin Ri.
158                for (int i = 0; i < k; ++i) {
159                        a[i] = (i*n)/k;
160                }
161
162                // (B)
163                int l, x;
164                for (int c = 0; c < k; ++c) {
165                        do {
166                                // Choose random x;
167                                x = 1 + random.nextInt(n);
168
169                                // determine range Rl;
170                                l = (x*k - 1)/n;
171                        } while (a[l] >= x); // accept or reject.
172
173                        ++a[l];
174                }
175                int s = k;
176
177                // (C) Move a[i] of nonempty bins to the left.
178                int m = 0, p = 0;
179                for (int i = 0; i < k; ++i) {
180                        if (a[i] == (i*n)/k) {
181                                a[i] = 0;
182                        } else {
183                                ++p;
184                                m = a[i];
185                                a[i] = 0;
186                                a[p - 1] = m;
187                        }
188                }
189
190                // (D) Determine l, set up space for Bl.
191                for (; p > 0; --p) {
192                        l = 1 + (a[p - 1]*k - 1)/n;
193                        final int ds = a[p - 1] - ((l - 1)*n)/k;
194                        a[p - 1] = 0;
195                        a[s - 1] = l;
196                        s -= ds;
197                }
198
199                // (E) If a[l] != 0, a new bin is to be processed.
200                int r = 0, m0 = 0;
201                for (int ll = 1; ll <= k; ++ll) {
202                        l = k + 1 - ll;
203
204                        if (a[l - 1] != 0) {
205                                r = l;
206                                m0 = 1 + ((a[l - 1] - 1)*n)/k;
207                                m = (a[l - 1]*n)/k - m0 + 1;
208                        }
209
210                        // (F) Choose a random x.
211                        x = m0 + random.nextInt(m);
212                        int i = l + 1;
213
214                        // (G) Check x against previously entered elements in bin;
215                        //     increment x as it jumps over elements <= x.
216                        while (i <= r && x >= a[i - 1]) {
217                                ++x;
218                                a[i - 2] = a[i - 1];
219                                ++i;
220                        }
221
222                        a[i - 2] = x;
223                        --m;
224                }
225
226                // Convert to zero based indexed arrays.
227                for (int i = 0; i < k; ++i) {
228                        --a[i];
229                }
230        }
231
232
233        private static void checkSubSet(final int n, final int k) {
234                if (k <= 0) {
235                        throw new IllegalArgumentException(format(
236                                "Subset size smaller or equal zero: %s", k
237                        ));
238                }
239                if (n < k) {
240                        throw new IllegalArgumentException(format(
241                                "n smaller than k: %s < %s.", n, k
242                        ));
243                }
244                if (!isMultiplicationSave(n, k)) {
245                        throw new IllegalArgumentException(format(
246                                "n*sub.length > Integer.MAX_VALUE (%s*%s = %s > %s)",
247                                n, k, (long)n*(long)k, Integer.MAX_VALUE
248                        ));
249                }
250        }
251
252}