From c8e3d899fa94638ab90fb8591a1b0a68c9077277 Mon Sep 17 00:00:00 2001 From: Jonathan Zheng Date: Thu, 16 Jan 2020 14:29:13 +0000 Subject: [PATCH] randomkit is good but not worth. Other optimizations needed --- lib/neatogen/CMakeLists.txt | 2 + lib/neatogen/Makefile.am | 4 +- lib/neatogen/randomkit.c | 245 ++++++++++++++++++++++++++++++++++++ lib/neatogen/randomkit.h | 150 ++++++++++++++++++++++ lib/neatogen/sgd.c | 11 +- 5 files changed, 409 insertions(+), 3 deletions(-) create mode 100644 lib/neatogen/randomkit.c create mode 100644 lib/neatogen/randomkit.h diff --git a/lib/neatogen/CMakeLists.txt b/lib/neatogen/CMakeLists.txt index ffdcf6c68..7e1780d20 100644 --- a/lib/neatogen/CMakeLists.txt +++ b/lib/neatogen/CMakeLists.txt @@ -46,6 +46,7 @@ set(sources stress.h voronoi.h sgd.h + randomkit.h # Source files adjust.c @@ -87,6 +88,7 @@ set(sources stress.c voronoi.c sgd.c + randomkit.c ) if (with_ipsepcola) diff --git a/lib/neatogen/Makefile.am b/lib/neatogen/Makefile.am index feea48a6a..487f43236 100644 --- a/lib/neatogen/Makefile.am +++ b/lib/neatogen/Makefile.am @@ -22,7 +22,7 @@ noinst_HEADERS = adjust.h edges.h geometry.h heap.h hedges.h info.h mem.h \ matrix_ops.h pca.h stress.h quad_prog_solver.h digcola.h \ overlap.h call_tri.h \ quad_prog_vpsc.h delaunay.h sparsegraph.h multispline.h fPQ.h \ - sgd.h + sgd.h randomkit.h IPSEPCOLA_SOURCES = constrained_majorization_ipsep.c \ mosek_quad_solve.c mosek_quad_solve.h quad_prog_vpsc.c @@ -39,6 +39,6 @@ libneatogen_C_la_SOURCES = adjust.c circuit.c edges.c geometry.c \ smart_ini_x.c constrained_majorization.c opt_arrangement.c \ overlap.c call_tri.c \ compute_hierarchy.c delaunay.c multispline.c $(WITH_IPSEPCOLA_SOURCES) \ - sgd.c + sgd.c randomkit.c EXTRA_DIST = $(IPSEPCOLA_SOURCES) gvneatogen.vcxproj* diff --git a/lib/neatogen/randomkit.c b/lib/neatogen/randomkit.c new file mode 100644 index 000000000..ba143a34b --- /dev/null +++ b/lib/neatogen/randomkit.c @@ -0,0 +1,245 @@ +/* Random kit 1.3 */ + +/* + * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org) + * + * The rk_random and rk_seed functions algorithms and the original design of + * the Mersenne Twister RNG: + * + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * 3. The names of its contributors may not be used to endorse or promote + * products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Original algorithm for the implementation of rk_interval function from + * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by + * Magnus Jonsson. + * + * Constants used in the rk_double implementation by Isaku Wada. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* static char const rcsid[] = + "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */ +#include +#include +#include +#include +#include +#include + +#ifdef _WIN32 +/* + * Windows + * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to + * detect the compiler used in distutils itself + */ +#if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND)) + +/* + * FIXME: ideally, we should set this to the real version of MSVCRT. We need + * something higher than 0x601 to enable _ftime64 and co + */ +#define __MSVCRT_VERSION__ 0x0700 +#include +#include + +/* + * mingw msvcr lib import wrongly export _ftime, which does not exist in the + * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which + * is available in those versions of the runtime + */ +#define _FTIME(x) _ftime64((x)) +#else +#include +#include +#define _FTIME(x) _ftime((x)) +#endif + +#ifndef RK_NO_WINCRYPT +/* Windows crypto */ +#ifndef _WIN32_WINNT +#define _WIN32_WINNT 0x0400 +#endif +#include +#include +#endif + +#else +/* Unix */ +#include +#include +#include +#endif + +#include "randomkit.h" + +#ifndef RK_DEV_URANDOM +#define RK_DEV_URANDOM "/dev/urandom" +#endif + +#ifndef RK_DEV_RANDOM +#define RK_DEV_RANDOM "/dev/random" +#endif + +char *rk_strerror[RK_ERR_MAX] = +{ + "no error", + "random device unvavailable" +}; + +void +rk_seed(unsigned long seed, rk_state *state) +{ + int pos; + seed &= 0xffffffffUL; + + /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ + for (pos = 0; pos < RK_STATE_LEN; pos++) { + state->key[pos] = seed; + seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL; + } + state->pos = RK_STATE_LEN; + state->gauss = 0; + state->has_gauss = 0; + state->has_binomial = 0; +} + +/* Magic Mersenne Twister constants */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL +#define UPPER_MASK 0x80000000UL +#define LOWER_MASK 0x7fffffffUL + +/* Slightly optimised reference implementation of the Mersenne Twister */ +unsigned long +rk_random(rk_state *state) +{ + unsigned long y; + + if (state->pos == RK_STATE_LEN) { + int i; + + for (i = 0; i < N - M; i++) { + y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + for (; i < N - 1; i++) { + y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); + state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A); + + state->pos = 0; + } + y = state->key[state->pos++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +long +rk_long(rk_state *state) +{ + return rk_ulong(state) >> 1; +} + +unsigned long +rk_ulong(rk_state *state) +{ +#if ULONG_MAX <= 0xffffffffUL + return rk_random(state); +#else + return (rk_random(state) << 32) | (rk_random(state)); +#endif +} + +unsigned long +rk_interval(unsigned long max, rk_state *state) +{ + unsigned long mask = max, value; + + if (max == 0) { + return 0; + } + /* Smallest bit mask >= max */ + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; +#if ULONG_MAX > 0xffffffffUL + mask |= mask >> 32; +#endif + + /* Search a random value in [0..mask] <= max */ +#if ULONG_MAX > 0xffffffffUL + if (max <= 0xffffffffUL) { + while ((value = (rk_random(state) & mask)) > max); + } + else { + while ((value = (rk_ulong(state) & mask)) > max); + } +#else + while ((value = (rk_ulong(state) & mask)) > max); +#endif + return value; +} + +double +rk_double(rk_state *state) +{ + /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ + long a = rk_random(state) >> 5, b = rk_random(state) >> 6; + return (a * 67108864.0 + b) / 9007199254740992.0; +} diff --git a/lib/neatogen/randomkit.h b/lib/neatogen/randomkit.h new file mode 100644 index 000000000..9676fa641 --- /dev/null +++ b/lib/neatogen/randomkit.h @@ -0,0 +1,150 @@ +/* Random kit 1.3 */ + +/* + * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: randomkit.h,v 1.24 2005/07/21 22:14:09 js Exp $ */ + +/* + * Typical use: + * + * { + * rk_state state; + * unsigned long seed = 1, random_value; + * + * rk_seed(seed, &state); // Initialize the RNG + * ... + * random_value = rk_random(&state); // Generate random values in [0..RK_MAX] + * } + * + * Instead of rk_seed, you can use rk_randomseed which will get a random seed + * from /dev/urandom (or the clock, if /dev/urandom is unavailable): + * + * { + * rk_state state; + * unsigned long random_value; + * + * rk_randomseed(&state); // Initialize the RNG with a random seed + * ... + * random_value = rk_random(&state); // Generate random values in [0..RK_MAX] + * } + */ + +/* + * Useful macro: + * RK_DEV_RANDOM: the device used for random seeding. + * defaults to "/dev/urandom" + */ + +#include + +#ifndef _RANDOMKIT_ +#define _RANDOMKIT_ + +#define RK_STATE_LEN 624 + +typedef struct rk_state_ +{ + unsigned long key[RK_STATE_LEN]; + int pos; + int has_gauss; /* !=0: gauss contains a gaussian deviate */ + double gauss; + + /* The rk_state structure has been extended to store the following + * information for the binomial generator. If the input values of n or p + * are different than nsave and psave, then the other parameters will be + * recomputed. RTK 2005-09-02 */ + + int has_binomial; /* !=0: following parameters initialized for + binomial */ + double psave; + long nsave; + double r; + double q; + double fm; + long m; + double p1; + double xm; + double xl; + double xr; + double c; + double laml; + double lamr; + double p2; + double p3; + double p4; + +} +rk_state; + +typedef enum { + RK_NOERR = 0, /* no error */ + RK_ENODEV = 1, /* no RK_DEV_RANDOM device */ + RK_ERR_MAX = 2 +} rk_error; + +/* error strings */ +extern char *rk_strerror[RK_ERR_MAX]; + +/* Maximum generated random value */ +#define RK_MAX 0xFFFFFFFFUL + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * Initialize the RNG state using the given seed. + */ +extern void rk_seed(unsigned long seed, rk_state *state); + +/* + * Returns a random unsigned long between 0 and RK_MAX inclusive + */ +extern unsigned long rk_random(rk_state *state); + +/* + * Returns a random long between 0 and LONG_MAX inclusive + */ +extern long rk_long(rk_state *state); + +/* + * Returns a random unsigned long between 0 and ULONG_MAX inclusive + */ +extern unsigned long rk_ulong(rk_state *state); + +/* + * Returns a random unsigned long between 0 and max inclusive. + */ +extern unsigned long rk_interval(unsigned long max, rk_state *state); + +/* + * Returns a random double between 0.0 and 1.0, 1.0 excluded. + */ +extern double rk_double(rk_state *state); + +#ifdef __cplusplus +} +#endif + +#endif /* _RANDOMKIT_ */ diff --git a/lib/neatogen/sgd.c b/lib/neatogen/sgd.c index ed63756bc..f7a69be7d 100644 --- a/lib/neatogen/sgd.c +++ b/lib/neatogen/sgd.c @@ -3,10 +3,13 @@ #include #include +#include "randomkit.h" + typedef struct term { node_t *i, *j; float d, w; } term; +// TODO: convert to just indices rather than node_t, because working with the auxiliary graph_t is slow float calculate_stress(term *terms, int n_terms) { float stress = 0; @@ -19,11 +22,16 @@ float calculate_stress(term *terms, int n_terms) { } return stress; } +// it is much faster to shuffle term rather than pointers to term, even though the swap is more expensive +static rk_state rstate; void fisheryates_shuffle(term *terms, int n_terms) { int i; for (i=n_terms-1; i>=1; i--) { // srand48() is called in neatoinit.c, so no need to seed here - int j = (int)(drand48() * (i+1)); + //int j = (int)(drand48() * (i+1)); + // TODO: better RNG because shuffling is eating up at least 50% of computation + + int j = rk_interval(i, &rstate); term temp = terms[i]; terms[i] = terms[j]; @@ -140,6 +148,7 @@ void sgd(graph_t *G, /* input graph */ start_timer(); } int t; + rk_seed(0, &rstate); for (t=0; t