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}