base.java
001 /*
002  * Java Genetic Algorithm Library (jenetics-3.5.0).
003  * Copyright (c) 2007-2016 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@gmx.at)
019  */
020 package org.jenetics.internal.math;
021 
022 import static java.lang.Double.doubleToLongBits;
023 import static java.lang.String.format;
024 import static java.util.Objects.requireNonNull;
025 
026 import java.util.Random;
027 
028 import org.jenetics.internal.util.require;
029 
030 import org.jenetics.util.RandomRegistry;
031 
032 /**
033  * This object contains mathematical helper functions.
034  *
035  * @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
036  * @since 1.0
037  * @version 3.0
038  */
039 public final class base {
040     private base() {require.noInstance();}
041 
042     /**
043      * <i>Clamping</i> a value between a pair of boundary values.
044      * <i>Note: using clamp with floating point numbers may give unexpected
045      * results if one of the values is {@code NaN}.</i>
046      *
047      * @param v the value to <i>clamp</i>
048      * @param lo the lower bound.
049      * @param hi the upper bound.
050      * @return The clamped value:
051      *        <ul>
052      *            <li>{@code lo if v < lo}</li>
053      *            <li>{@code hi if hi < v}</li>
054      *            <li>{@code otherwise, v}</li>
055      *        </ul>
056      */
057     public static double clamp(final double v, final double lo, final double hi) {
058         return v < lo ? lo : (v > hi ? hi : v);
059     }
060 
061     /**
062      * Return the <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place">ULP</a>
063      * distance of the given two double values.
064      *
065      * @param a first double.
066      * @param b second double.
067      * @return the ULP distance.
068      * @throws ArithmeticException if the distance doesn't fit in a long value.
069      */
070     public static long ulpDistance(final double a, final double b) {
071         return Math.subtractExact(ulpPosition(a), ulpPosition(b));
072     }
073 
074     /**
075      * Calculating the <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place">ULP</a>
076      * position of a double number.
077      *
078      * <pre>{@code
079      * double a = 0.0;
080      * for (int i = 0; i < 10; ++i) {
081      *     a = Math.nextAfter(a, Double.POSITIVE_INFINITY);
082      * }
083      *
084      * for (int i = 0; i < 19; ++i) {
085      *     a = Math.nextAfter(a, Double.NEGATIVE_INFINITY);
086      *     System.out.println(
087      *          a + "\t" + ulpPosition(a) + "\t" + ulpDistance(0.0, a)
088      *     );
089      * }
090      * }</pre>
091      *
092      * The code fragment above will create the following output:
093      * <pre>
094      *   4.4E-323    9  9
095      *   4.0E-323    8  8
096      *   3.5E-323    7  7
097      *   3.0E-323    6  6
098      *   2.5E-323    5  5
099      *   2.0E-323    4  4
100      *   1.5E-323    3  3
101      *   1.0E-323    2  2
102      *   4.9E-324    1  1
103      *   0.0         0  0
104      *  -4.9E-324   -1  1
105      *  -1.0E-323   -2  2
106      *  -1.5E-323   -3  3
107      *  -2.0E-323   -4  4
108      *  -2.5E-323   -5  5
109      *  -3.0E-323   -6  6
110      *  -3.5E-323   -7  7
111      *  -4.0E-323   -8  8
112      *  -4.4E-323   -9  9
113      * </pre>
114      *
115      * @param a the double number.
116      * @return the ULP position.
117      */
118     public static long ulpPosition(final double a) {
119         long t = doubleToLongBits(a);
120         if (t < 0) {
121             t = Long.MIN_VALUE - t;
122         }
123         return t;
124     }
125 
126     /**
127      * Selects a random subset of size {@code k} from a set of size {@code n}.
128      *
129      * @see #subset(int, int[])
130      *
131      * @param n the size of the set.
132      * @param k the size of the subset.
133      * @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
134      *          {@code n*k} will cause an integer overflow.
135      * @return the subset array.
136      */
137     public static int[] subset(final int n, final int k) {
138         return subset(n, k, RandomRegistry.getRandom());
139     }
140 
141     /**
142      * Selects a random subset of size {@code k} from a set of size {@code n}.
143      *
144      * @see #subset(int, int[], Random)
145      *
146      * @param n the size of the set.
147      * @param k the size of the subset.
148      * @param random the random number generator used.
149      * @throws NullPointerException if {@code random} is {@code null}.
150      * @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
151      *         {@code n*k} will cause an integer overflow.
152      * @return the subset array.
153      */
154     public static int[] subset(final int n, final int k, final Random random) {
155         requireNonNull(random, "Random");
156         if (k <= 0) {
157             throw new IllegalArgumentException(format(
158                     "Subset size smaller or equal zero: %s", k
159                 ));
160         }
161         if (n < k) {
162             throw new IllegalArgumentException(format(
163                     "n smaller than k: %s < %s.", n, k
164                 ));
165         }
166 
167         final int[] sub = new int[k];
168         subset(n, sub,random);
169         return sub;
170     }
171 
172     /**
173      * <p>
174      * Selects a random subset of size {@code sub.length} from a set of size
175      * {@code n}.
176      * </p>
177      *
178      * <p>
179      * <em>Authors:</em>
180      *      FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. This
181      *      version based on the  C++ version by John Burkardt.
182      * </p>
183      *
184      * <p><em><a href="https://people.scs.fsu.edu/~burkardt/c_src/subset/subset.html">
185      *  Reference:</a></em>
186      *      Albert Nijenhuis, Herbert Wilf,
187      *      Combinatorial Algorithms for Computers and Calculators,
188      *      Second Edition,
189      *      Academic Press, 1978,
190      *      ISBN: 0-12-519260-6,
191      *      LC: QA164.N54.
192      * </p>
193      *
194      * @param n the size of the set.
195      * @param sub the sub set array.
196      * @throws NullPointerException if {@code sub} is {@code null}.
197      * @throws IllegalArgumentException if {@code n < sub.length},
198      *         {@code sub.length == 0} or {@code n*sub.length} will cause an
199      *         integer overflow.
200      */
201     public static void subset(final int n, final int sub[]) {
202         subset(n, sub, RandomRegistry.getRandom());
203     }
204 
205     /**
206      * <p>
207      * Selects a random subset of size {@code sub.length} from a set of size
208      * {@code n}.
209      * </p>
210      *
211      * <p>
212      * <em>Authors:</em>
213      *      FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. This
214      *      version based on the  C++ version by John Burkardt.
215      * </p>
216      *
217      * <p><em><a href="https://people.scs.fsu.edu/~burkardt/c_src/subset/subset.html">
218      *  Reference:</a></em>
219      *      Albert Nijenhuis, Herbert Wilf,
220      *      Combinatorial Algorithms for Computers and Calculators,
221      *      Second Edition,
222      *      Academic Press, 1978,
223      *      ISBN: 0-12-519260-6,
224      *      LC: QA164.N54.
225      * </p>
226      *
227      * @param n the size of the set.
228      * @param sub the sub set array.
229      * @param random the random number generator used.
230      * @return the sub-set array for the given parameter
231      * @throws NullPointerException if {@code sub} or {@code random} is
232      *         {@code null}.
233      * @throws IllegalArgumentException if {@code n < sub.length},
234      *         {@code sub.length == 0} or {@code n*sub.length} will cause an
235      *         integer overflow.
236      */
237     public static int[] subset(final int n, final int sub[], final Random random) {
238         requireNonNull(random, "Random");
239         requireNonNull(sub, "Sub set array");
240 
241         final int k = sub.length;
242         checkSubSet(n, k);
243 
244         if (sub.length == n) {
245             for (int i = 0; i < k; ++i) {
246                 sub[i] = i;
247             }
248             return sub;
249         }
250 
251         for (int i = 0; i < k; ++i) {
252             sub[i] = (i*n)/k;
253         }
254 
255         int l = 0;
256         int ix = 0;
257         for (int i = 0; i < k; ++i) {
258             do {
259                 ix = nextInt(random, 1, n);
260                 l = (ix*k - 1)/n;
261             } while (sub[l] >= ix);
262 
263             sub[l] = sub[l] + 1;
264         }
265 
266         int m = 0;
267         int ip = 0;
268         int is = k;
269         for (int i = 0; i < k; ++i) {
270             m = sub[i];
271             sub[i] = 0;
272 
273             if (m != (i*n)/k) {
274                 ip = ip + 1;
275                 sub[ip - 1] = m;
276             }
277         }
278 
279         int ihi = ip;
280         int ids = 0;
281         for (int i = 1; i <= ihi; ++i) {
282             ip = ihi + 1 - i;
283             l = 1 + (sub[ip - 1]*k - 1)/n;
284             ids = sub[ip - 1] - ((l - 1)*n)/k;
285             sub[ip - 1] = 0;
286             sub[is - 1] = l;
287             is = is - ids;
288         }
289 
290         int ir = 0;
291         int m0 = 0;
292         for (int ll = 1; ll <= k; ++ll) {
293             l = k + 1 - ll;
294 
295             if (sub[l - 1] != 0) {
296                 ir = l;
297                 m0 = 1 + ((sub[l - 1] - 1)*n)/k;
298                 m = (sub[l-1]*n)/k - m0 + 1;
299             }
300 
301             ix = nextInt(random, m0, m0 + m - 1);
302 
303             int i = l + 1;
304             while (i <= ir && ix >= sub[i - 1]) {
305                 ix = ix + 1;
306                 sub[ i- 2] = sub[i - 1];
307                 i = i + 1;
308             }
309 
310             sub[i - 2] = ix;
311             --m;
312         }
313 
314         return sub;
315     }
316 
317     public static void checkSubSet(final int n, final int k) {
318         if (k <= 0) {
319             throw new IllegalArgumentException(format(
320                 "Subset size smaller or equal zero: %s", k
321             ));
322         }
323         if (n < k) {
324             throw new IllegalArgumentException(format(
325                 "n smaller than k: %s < %s.", n, k
326             ));
327         }
328         if (!arithmetic.isMultiplicationSave(n, k)) {
329             throw new IllegalArgumentException(format(
330                 "n*sub.length > Integer.MAX_VALUE (%s*%s = %s > %s)",
331                 n, k, (long)n*(long)k, Integer.MAX_VALUE
332             ));
333         }
334     }
335 
336     private static int nextInt(final Random random, final int a, final int b) {
337         return a == b ? a - 1 : random.nextInt(b - a) + a;
338     }
339 
340 }