diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 048e0a3b18..25f1af9c3b 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -448,22 +448,10 @@ set(SCOPMATH_FILES_LIST set(NRNMPI_FILES_LIST nrnmpi.cpp bbsmpipack.cpp mpispike.cpp) set(NRNGNU_FILES_LIST - ACG.cpp - Binomial.cpp - DiscUnif.cpp - Erlang.cpp - Geom.cpp HypGeom.cpp - LogNorm.cpp - MLCG.cpp - NegExp.cpp - Normal.cpp - Poisson.cpp RNG.cpp - Random.cpp - RndInt.cpp - Uniform.cpp - Weibull.cpp) + RNG_random123.cpp + Random.cpp) # nrnpython sources (only if ${NRN_ENABLE_PYTHON_DYNAMIC} is OFF} set(NRNPYTHON_FILES_LIST diff --git a/src/gnu/ACG.cpp b/src/gnu/ACG.cpp deleted file mode 100755 index ecba08420b..0000000000 --- a/src/gnu/ACG.cpp +++ /dev/null @@ -1,293 +0,0 @@ -#include <../../nrnconf.h> -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1989 Free Software Foundation - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. -*/ - -#include -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -// -// This is an extension of the older implementation of Algorithm M -// which I previously supplied. The main difference between this -// version and the old code are: -// -// + Andres searched high & low for good constants for -// the LCG. -// -// + theres more bit chopping going on. -// -// The following contains his comments. -// -// agn@UNH.CS.CMU.EDU sez.. -// -// The generator below is based on 2 well known -// methods: Linear Congruential (LCGs) and Additive -// Congruential generators (ACGs). -// -// The LCG produces the longest possible sequence -// of 32 bit random numbers, each being unique in -// that sequence (it has only 32 bits of state). -// It suffers from 2 problems: a) Independence -// isnt great, that is the (n+1)th number is -// somewhat related to the preceding one, unlike -// flipping a coin where knowing the past outcomes -// dont help to predict the next result. b) -// Taking parts of a LCG generated number can be -// quite non-random: for example, looking at only -// the least significant byte gives a permuted -// 8-bit counter (that has a period length of only -// 256). The advantage of an LCA is that it is -// perfectly uniform when run for the entire period -// length (and very uniform for smaller sequences -// too, if the parameters are chosen carefully). -// -// ACGs have extremly long period lengths and -// provide good independence. Unfortunately, -// uniformity isnt not too great. Furthermore, I -// didnt find any theoretically analysis of ACGs -// that addresses uniformity. -// -// The RNG given below will return numbers -// generated by an LCA that are permuted under -// control of a ACG. 2 permutations take place: the -// 4 bytes of one LCG generated number are -// subjected to one of 16 permutations selected by -// 4 bits of the ACG. The permutation a such that -// byte of the result may come from each byte of -// the LCG number. This effectively destroys the -// structure within a word. Finally, the sequence -// of such numbers is permuted within a range of -// 256 numbers. This greatly improves independence. -// -// -// Algorithm M as describes in Knuths "Art of Computer Programming", -// Vol 2. 1969 -// is used with a linear congruential generator (to get a good uniform -// distribution) that is permuted with a Fibonacci additive congruential -// generator to get good independence. -// -// Bit, byte, and word distributions were extensively tested and pass -// Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity -// assumption holds with probability > 0.999) -// -// Run-up tests for on 7E8 numbers confirm independence with -// probability > 0.97. -// -// Plotting random points in 2d reveals no apparent structure. -// -// Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i), -// i=1..512) -// results in no obvious structure (A(i) ~ const). -// -// Except for speed and memory requirements, this generator outperforms -// random() for all tests. (random() scored rather low on uniformity tests, -// while independence test differences were less dramatic). -// -// AGN would like to.. -// thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help. -// -// And I would (DGC) would like to thank Donald Kunth for AGN for letting me -// use his extensions in this implementation. -// - -// -// Part of the table on page 28 of Knuth, vol II. This allows us -// to adjust the size of the table at the expense of shorter sequences. -// - -static int randomStateTable[][3] = { -{3,7,16}, {4,9, 32}, {3,10, 32}, {1,11, 32}, {1,15,64}, {3,17,128}, -{7,18,128}, {3,20,128}, {2,21, 128}, {1,22, 128}, {5,23, 128}, {3,25, 128}, -{2,29, 128}, {3,31, 128}, {13,33, 256}, {2,35, 256}, {11,36, 256}, -{14,39,256}, {3,41,256}, {9,49,256}, {3,52,256}, {24,55,256}, {7,57, 256}, -{19,58,256}, {38,89,512}, {17,95,512}, {6,97,512}, {11,98,512}, {-1,-1,-1} }; - -// -// spatial permutation table -// RANDOM_PERM_SIZE must be a power of two -// - -#define RANDOM_PERM_SIZE 64 -uint32_t randomPermutations[RANDOM_PERM_SIZE] = { -0xffffffff, 0x00000000, 0x00000000, 0x00000000, // 3210 -0x0000ffff, 0x00ff0000, 0x00000000, 0xff000000, // 2310 -0xff0000ff, 0x0000ff00, 0x00000000, 0x00ff0000, // 3120 -0x00ff00ff, 0x00000000, 0xff00ff00, 0x00000000, // 1230 - -0xffff0000, 0x000000ff, 0x00000000, 0x0000ff00, // 3201 -0x00000000, 0x00ff00ff, 0x00000000, 0xff00ff00, // 2301 -0xff000000, 0x00000000, 0x000000ff, 0x00ffff00, // 3102 -0x00000000, 0x00000000, 0x00000000, 0xffffffff, // 2103 - -0xff00ff00, 0x00000000, 0x00ff00ff, 0x00000000, // 3012 -0x0000ff00, 0x00000000, 0x00ff0000, 0xff0000ff, // 2013 -0x00000000, 0x00000000, 0xffffffff, 0x00000000, // 1032 -0x00000000, 0x0000ff00, 0xffff0000, 0x000000ff, // 1023 - -0x00000000, 0xffffffff, 0x00000000, 0x00000000, // 0321 -0x00ffff00, 0xff000000, 0x00000000, 0x000000ff, // 0213 -0x00000000, 0xff000000, 0x0000ffff, 0x00ff0000, // 0132 -0x00000000, 0xff00ff00, 0x00000000, 0x00ff00ff // 0123 -}; - -// -// SEED_TABLE_SIZE must be a power of 2 -// -#define SEED_TABLE_SIZE 32 -static uint32_t seedTable[SEED_TABLE_SIZE] = { -0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b, -0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf, -0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706, -0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10, -0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a, -0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32, -0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f, -0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf -}; - -// -// The LCG used to scramble the ACG -// -// -// LC-parameter selection follows recommendations in -// "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi. -// -// LC_A = 251^2, ~= sqrt(2^32) = 66049 -// LC_C = result of a long trial & error series = 3907864577 -// - -static const uint32_t LC_A = 66049; -static const uint32_t LC_C = 3907864577U; -static inline uint32_t LCG(uint32_t x) -{ - return( x * LC_A + LC_C ); -} - - -ACG::ACG(uint32_t seed, int size) -{ - int l; - initialSeed = seed; - // - // Determine the size of the state table - // - - for (l = 0; - randomStateTable[l][0] != -1 && randomStateTable[l][1] < size; - l++); - - if (randomStateTable[l][1] == -1) { - l--; - } - - initialTableEntry = l; - - stateSize = randomStateTable[ initialTableEntry ][ 1 ]; - auxSize = randomStateTable[ initialTableEntry ][ 2 ]; - - // - // Allocate the state table & the auxillary table in a single malloc - // - - state = new uint32_t[stateSize + auxSize]; - auxState = &state[stateSize]; - - reset(); -} - -// -// Initialize the state -// -void -ACG::reset() -{ - uint32_t u; - - if (initialSeed < SEED_TABLE_SIZE) { - u = seedTable[ initialSeed ]; - } else { - u = initialSeed ^ seedTable[ initialSeed & (SEED_TABLE_SIZE-1) ]; - } - - - j = randomStateTable[ initialTableEntry ][ 0 ] - 1; - k = randomStateTable[ initialTableEntry ][ 1 ] - 1; - - int i; - for(i = 0; i < stateSize; i++) { - state[i] = u = LCG(u); - } - - for (i = 0; i < auxSize; i++) { - auxState[i] = u = LCG(u); - } - - k = u % stateSize; - int tailBehind = (stateSize - randomStateTable[ initialTableEntry ][ 0 ]); - j = k - tailBehind; - if (j < 0) { - j += stateSize; - } - - lcgRecurr = u; - - assert(sizeof(double) == 2 * sizeof(int32_t)); -} - -ACG::~ACG() -{ - if (state) delete [] state; - state = 0; - // don't delete auxState, it's really an alias for state. -} - -// -// Returns 32 bits of random information. -// - -uint32_t -ACG::asLong() -{ - uint32_t result = state[k] + state[j]; - state[k] = result; - j = (j <= 0) ? (stateSize-1) : (j-1); - k = (k <= 0) ? (stateSize-1) : (k-1); - - short int auxIndex = (result >> 24) & (auxSize - 1); - uint32_t auxACG = auxState[auxIndex]; - auxState[auxIndex] = lcgRecurr = LCG(lcgRecurr); - - // - // 3c is a magic number. We are doing four masks here, so we - // do not want to run off the end of the permutation table. - // This insures that we have always got four entries left. - // - uint32_t *perm = & randomPermutations[result & 0x3c]; - - result = *(perm++) & auxACG; - result |= *(perm++) & ((auxACG << 24) - | ((auxACG >> 8)& 0xffffff)); - result |= *(perm++) & ((auxACG << 16) - | ((auxACG >> 16) & 0xffff)); - result |= *(perm++) & ((auxACG << 8) - | ((auxACG >> 24) & 0xff)); - - return(result); -} diff --git a/src/gnu/ACG.h b/src/gnu/ACG.h deleted file mode 100755 index 8cee68fafa..0000000000 --- a/src/gnu/ACG.h +++ /dev/null @@ -1,68 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _ACG_h -#define _ACG_h 1 - -#include -#include -#ifdef __GNUG__ -//#pragma interface -#endif - -// -// Additive number generator. This method is presented in Volume II -// of The Art of Computer Programming by Knuth. I've coded the algorithm -// and have added the extensions by Andres Nowatzyk of CMU to randomize -// the result of algorithm M a bit by using an LCG & a spatial -// permutation table. -// -// The version presented uses the same constants for the LCG that Andres -// uses (chosen by trial & error). The spatial permutation table is -// the same size (it's based on word size). This is for 32-bit words. -// -// The ``auxillary table'' used by the LCG table varies in size, and -// is chosen to be the the smallest power of two which is larger than -// twice the size of the state table. -// - -class ACG : public RNG { - - uint32_t initialSeed; // used to reset generator - int initialTableEntry; - - uint32_t *state; - uint32_t *auxState; - short stateSize; - short auxSize; - uint32_t lcgRecurr; - short j; - short k; - -protected: - -public: - ACG(uint32_t seed = 0, int size = 55); - virtual ~ACG(); - // - // Return a long-words word of random bits - // - virtual uint32_t asLong(); - virtual void reset(); -}; - -#endif diff --git a/src/gnu/Binomial.cpp b/src/gnu/Binomial.cpp deleted file mode 100755 index 3794adc75a..0000000000 --- a/src/gnu/Binomial.cpp +++ /dev/null @@ -1,34 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double Binomial::operator()() -{ - int s = 0; - for (int i = 0; i < pN; i++) { - if (pGenerator -> asDouble() < pU) { - s++; - } - } - return(double(s)); -} - diff --git a/src/gnu/Binomial.h b/src/gnu/Binomial.h deleted file mode 100755 index 49d9c9170f..0000000000 --- a/src/gnu/Binomial.h +++ /dev/null @@ -1,55 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Binomial_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Binomial_h 1 - -#include - -class Binomial: public Random { -protected: - int pN; - double pU; -public: - Binomial(int n, double u, RNG *gen); - - int n(); - int n(int xn); - - double u(); - double u(int xu); - - virtual double operator()(); - -}; - - -inline Binomial::Binomial(int n, double u, RNG *gen) -: Random(gen){ - pN = n; pU = u; -} - -inline int Binomial::n() { return pN; } -inline int Binomial::n(int xn) { int tmp = pN; pN = xn; return tmp; } - -inline double Binomial::u() { return pU; } -inline double Binomial::u(int xu) { double tmp = pU; pU = xu; return tmp; } - -#endif diff --git a/src/gnu/DiscUnif.cpp b/src/gnu/DiscUnif.cpp deleted file mode 100755 index 72a9df8222..0000000000 --- a/src/gnu/DiscUnif.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double DiscreteUniform::operator()() -{ - long tmp = long(floor(delta * pGenerator -> asDouble())); - return( double(pLow + tmp) ); -} - diff --git a/src/gnu/DiscUnif.h b/src/gnu/DiscUnif.h deleted file mode 100755 index 5fa98eb4a9..0000000000 --- a/src/gnu/DiscUnif.h +++ /dev/null @@ -1,72 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _DiscreteUniform_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _DiscreteUniform_h 1 - -#include - -// -// The interval [lo..hi) -// - -class DiscreteUniform: public Random { - long pLow; - long pHigh; - double delta; -public: - DiscreteUniform(long low, long high, RNG *gen); - - long low(); - long low(long x); - long high(); - long high(long x); - - virtual double operator()(); -}; - - -inline DiscreteUniform::DiscreteUniform(long low, long high, RNG *gen) -: Random(gen) -{ - pLow = (low < high) ? low : high; - pHigh = (low < high) ? high : low; - delta = (pHigh - pLow) + 1; -} - -inline long DiscreteUniform::low() { return pLow; } - -inline long DiscreteUniform::low(long x) { - long tmp = pLow; - pLow = x; - delta = (pHigh - pLow) + 1; - return tmp; -} - -inline long DiscreteUniform::high() { return pHigh; } - -inline long DiscreteUniform::high(long x) { - long tmp = pHigh; - pHigh = x; - delta = (pHigh - pLow) + 1; - return tmp; -} - -#endif diff --git a/src/gnu/Erlang.cpp b/src/gnu/Erlang.cpp deleted file mode 100755 index f9638df073..0000000000 --- a/src/gnu/Erlang.cpp +++ /dev/null @@ -1,32 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double Erlang::operator()() -{ - double prod = 1.0; - - for (int i = 0; i < k; i++) { - prod *= pGenerator -> asDouble(); - } - return(-log(prod)/a); -} diff --git a/src/gnu/Erlang.h b/src/gnu/Erlang.h deleted file mode 100755 index d6c8ccd5e3..0000000000 --- a/src/gnu/Erlang.h +++ /dev/null @@ -1,68 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Erlang_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Erlang_h 1 - -#include - -class Erlang: public Random { -protected: - double pMean; - double pVariance; - int k; - double a; - void setState(); -public: - Erlang(double mean, double variance, RNG *gen); - - double mean(); - double mean(double x); - double variance(); - double variance(double x); - - virtual double operator()(); - -}; - - -inline void Erlang::setState() { - k = int( (pMean * pMean ) / pVariance + 0.5 ); - k = (k > 0) ? k : 1; - a = k / pMean; -} - -inline Erlang::Erlang(double mean, double variance, RNG *gen) : Random(gen) -{ - pMean = mean; pVariance = variance; - setState(); -} - -inline double Erlang::mean() { return pMean; } -inline double Erlang::mean(double x) { - double tmp = pMean; pMean = x; setState(); return tmp; -}; - -inline double Erlang::variance() { return pVariance; } -inline double Erlang::variance(double x) { - double tmp = pVariance; pVariance = x; setState(); return tmp; -} - -#endif diff --git a/src/gnu/Geom.cpp b/src/gnu/Geom.cpp deleted file mode 100755 index d2883f5469..0000000000 --- a/src/gnu/Geom.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double Geometric::operator()() -{ - int samples; - for (samples = 1; pGenerator -> asDouble() < pMean; samples++); - return((double) samples); -} - diff --git a/src/gnu/Geom.h b/src/gnu/Geom.h deleted file mode 100755 index 231744fe41..0000000000 --- a/src/gnu/Geom.h +++ /dev/null @@ -1,52 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Geometric_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Geometric_h - -#include - -class Geometric: public Random { -protected: - double pMean; -public: - Geometric(double mean, RNG *gen); - - double mean(); - double mean(double x); - - virtual double operator()(); - -}; - - -inline Geometric::Geometric(double mean, RNG *gen) : Random(gen) -{ - pMean = mean; -} - - -inline double Geometric::mean() { return pMean; } -inline double Geometric::mean(double x) { - double tmp = pMean; pMean = x; return tmp; -} - - -#endif diff --git a/src/gnu/HypGeom.cpp b/src/gnu/HypGeom.cpp index 9dc1ba5d2b..ef23c198b6 100755 --- a/src/gnu/HypGeom.cpp +++ b/src/gnu/HypGeom.cpp @@ -24,7 +24,7 @@ Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. double HyperGeometric::operator()() { - double d = (pGenerator -> asDouble() > pP) ? (1.0 - pP) : (pP); - return(-pMean * log(pGenerator -> asDouble()) / (2.0 * d) ); + double d = ((*(generator()->get_random_generator()))() > pP) ? (1.0 - pP) : (pP); + return(-pMean * log((*(generator()->get_random_generator()))()) / (2.0 * d) ); } diff --git a/src/gnu/HypGeom.h b/src/gnu/HypGeom.h index ea99afe70b..6c585857f0 100755 --- a/src/gnu/HypGeom.h +++ b/src/gnu/HypGeom.h @@ -23,7 +23,7 @@ Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. #include -class HyperGeometric: public Random { +class HyperGeometric: public Random_random123 { protected: double pMean; double pVariance; @@ -31,7 +31,7 @@ class HyperGeometric: public Random { void setState(); public: - HyperGeometric(double mean, double variance, RNG *gen); + HyperGeometric(double mean, double variance, RNG_random123 *gen); double mean(); double mean(double x); @@ -47,8 +47,8 @@ inline void HyperGeometric::setState() { pP = 0.5 * (1.0 - sqrt((z - 1.0) / ( z + 1.0 ))); } -inline HyperGeometric::HyperGeometric(double mean, double variance, RNG *gen) -: Random(gen) { +inline HyperGeometric::HyperGeometric(double mean, double variance, RNG_random123 *gen) +: Random_random123(gen) { pMean = mean; pVariance = variance; setState(); } diff --git a/src/gnu/LogNorm.cpp b/src/gnu/LogNorm.cpp deleted file mode 100755 index 14f75f3fca..0000000000 --- a/src/gnu/LogNorm.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -#include - -#ifndef M_E -#define M_E 2.71828182845904523536 -#endif - -// -// See Simulation, Modelling & Analysis by Law & Kelton, pp260 -// -// - -double LogNormal::operator()() -{ - return( pow(M_E, this->Normal::operator()() ) ); -} - - diff --git a/src/gnu/LogNorm.h b/src/gnu/LogNorm.h deleted file mode 100755 index ccc94ce4be..0000000000 --- a/src/gnu/LogNorm.h +++ /dev/null @@ -1,78 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _LogNormal_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _LogNormal_h - -#include - -class LogNormal: public Normal { -protected: - double logMean; - double logVariance; - void setState(); -public: - LogNormal(double mean, double variance, RNG *gen); - double mean(); - double mean(double x); - double variance(); - double variance(double x); - virtual double operator()(); -}; - - -inline void LogNormal::setState() -{ - double m2 = logMean * logMean; - pMean = log(m2 / sqrt(logVariance + m2) ); -// from ch@heike.informatik.uni-dortmund.de: -// (was pVariance = log((sqrt(logVariance + m2)/m2 )); ) - pStdDev = sqrt(log((logVariance + m2)/m2 )); -} - -inline LogNormal::LogNormal(double mean, double variance, RNG *gen) - : Normal(mean, variance, gen) -{ - logMean = mean; - logVariance = variance; - setState(); -} - -inline double LogNormal::mean() { - return logMean; -} - -inline double LogNormal::mean(double x) -{ - double t=logMean; logMean = x; setState(); - return t; -} - -inline double LogNormal::variance() { - return logVariance; -} - -inline double LogNormal::variance(double x) -{ - double t=logVariance; logVariance = x; setState(); - return t; -} - -#endif diff --git a/src/gnu/MLCG.cpp b/src/gnu/MLCG.cpp deleted file mode 100755 index aa8b10f4b8..0000000000 --- a/src/gnu/MLCG.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#include <../../nrnconf.h> -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1989 Free Software Foundation - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -// -// SEED_TABLE_SIZE must be a power of 2 -// - - -#define SEED_TABLE_SIZE 32 - -static int32_t seedTable[SEED_TABLE_SIZE] = { - /* 0xbdcc47e5 */ -1110685723, /* 0x54aea45d */ 1420731485, - /* 0xec0df859 */ -334628775, /* 0xda84637b */ -628857989, - /* 0xc8c6cb4f */ -926495921, /* 0x35574b01 */ 894913281, - /* 0x28260b7d */ 673581949, /* 0x0d07fdbf */ 218627519, - /* 0x9faaeeb0 */ -1616187728, /* 0x613dd169 */ 1631441257, - /* 0x5ce2d818 */ 1558370328, /* 0x85b9e706 */ -2051414266, - /* 0xab2469db */ -1423676965, /* 0xda02b0dc */ -637357860, - /* 0x45c60d6e */ 1170607470, /* 0xffe49d10 */ -1794800, - /* 0x7224fea3 */ 1915027107, /* 0xf9684fc9 */ -110604343, - /* 0xfc7ee074 */ -58793868, /* 0x326ce92a */ 845998378, - /* 0x366d13b5 */ 913118133, /* 0x17aaa731 */ 397059889, - /* 0xeb83a675 */ -343693707, /* 0x7781cb32 */ 2004994866, - /* 0x4ec7c92d */ 1321716013, /* 0x7f187521 */ 2132309281, - /* 0x2cf346b4 */ 754140852, /* 0xad13310f */ -1391251185, - /* 0xb89cff2b */ -1197670613, /* 0x12164de1 */ 303451617, - /* 0xa865168d */ -1469770099, /* 0x32b56cdf */ 850750687 -}; - -MLCG::MLCG(int32_t seed1, int32_t seed2) -{ - initialSeedOne = seed1; - initialSeedTwo = seed2; - reset(); -} - -void -MLCG::reset() -{ - int32_t seed1 = initialSeedOne; - int32_t seed2 = initialSeedTwo; - - // Most people pick stupid seed numbers that do not have enough - // bits. In this case, if they pick a small seed number, we - // map that to a specific seed. - - if (seed1 < 0) { - seed1 = (seed1 + 2147483561); - seed1 = (seed1 < 0) ? -seed1 : seed1; - } - - if (seed2 < 0) { - seed2 = (seed2 + 2147483561); - seed2 = (seed2 < 0) ? -seed2 : seed2; - } - - if (seed1 > -1 && seed1 < SEED_TABLE_SIZE) { - seedOne = seedTable[seed1]; - } else { - seedOne = seed1 ^ seedTable[seed1 & (SEED_TABLE_SIZE-1)]; - } - - if (seed2 > -1 && seed2 < SEED_TABLE_SIZE) { - seedTwo = seedTable[seed2]; - } else { - seedTwo = seed2 ^ seedTable[ seed2 & (SEED_TABLE_SIZE-1) ]; - } - seedOne = (seedOne % 2147483561) + 1; - seedTwo = (seedTwo % 2147483397) + 1; -} - -uint32_t MLCG::asLong() -{ - int32_t k = seedOne % 53668; - - seedOne = 40014 * (seedOne-k * 53668) - k * 12211; - if (seedOne < 0) { - seedOne += 2147483563; - } - - k = seedTwo % 52774; - seedTwo = 40692 * (seedTwo - k * 52774) - k * 3791; - if (seedTwo < 0) { - seedTwo += 2147483399; - } - - int32_t z = seedOne - seedTwo; - if (z < 1) { - z += 2147483562; - } - return( (unsigned long) z); -} - diff --git a/src/gnu/MLCG.h b/src/gnu/MLCG.h deleted file mode 100755 index e0ff5c4b62..0000000000 --- a/src/gnu/MLCG.h +++ /dev/null @@ -1,87 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _MLCG_h -#define _MLCG_h 1 -#ifdef __GNUG__ -//#pragma interface -#endif - -#include -#include - -// -// Multiplicative Linear Conguential Generator -// - -class MLCG : public RNG { - int32_t initialSeedOne; - int32_t initialSeedTwo; - int32_t seedOne; - int32_t seedTwo; - -protected: - -public: - MLCG(int32_t seed1 = 0, int32_t seed2 = 1); - // - // Return a long-words word of random bits - // - virtual uint32_t asLong(); - virtual void reset(); - int32_t seed1(); - void seed1(int32_t); - int32_t seed2(); - void seed2(int32_t); - void reseed(int32_t, int32_t); -}; - -inline int32_t -MLCG::seed1() -{ - return(seedOne); -} - -inline void -MLCG::seed1(int32_t s) -{ - initialSeedOne = s; - reset(); -} - -inline int32_t -MLCG::seed2() -{ - return(seedTwo); -} - -inline void -MLCG::seed2(int32_t s) -{ - initialSeedTwo = s; - reset(); -} - -inline void -MLCG::reseed(int32_t s1, int32_t s2) -{ - initialSeedOne = s1; - initialSeedTwo = s2; - reset(); -} - -#endif diff --git a/src/gnu/NegExp.cpp b/src/gnu/NegExp.cpp deleted file mode 100755 index caf5afc913..0000000000 --- a/src/gnu/NegExp.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double NegativeExpntl::operator()() -{ - return(-pMean * log(pGenerator -> asDouble())); -} - diff --git a/src/gnu/NegExp.h b/src/gnu/NegExp.h deleted file mode 100755 index 1be7d68e8b..0000000000 --- a/src/gnu/NegExp.h +++ /dev/null @@ -1,55 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _NegativeExpntl_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _NegativeExpntl_h 1 - - -// -// Negative Exponential Random Numbers -// -// - -#include - -class NegativeExpntl: public Random { -protected: - double pMean; -public: - NegativeExpntl(double xmean, RNG *gen); - double mean(); - double mean(double x); - - virtual double operator()(); -}; - - -inline NegativeExpntl::NegativeExpntl(double xmean, RNG *gen) -: Random(gen) { - pMean = xmean; -} - -inline double NegativeExpntl::mean() { return pMean; } -inline double NegativeExpntl::mean(double x) { - double t = pMean; pMean = x; - return t; -} - -#endif diff --git a/src/gnu/Normal.cpp b/src/gnu/Normal.cpp deleted file mode 100755 index 25aab78752..0000000000 --- a/src/gnu/Normal.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include -// -// See Simulation, Modelling & Analysis by Law & Kelton, pp259 -// -// This is the ``polar'' method. -// - -double Normal::operator()() -{ - - if (haveCachedNormal == 1) { - haveCachedNormal = 0; - return(cachedNormal * pStdDev + pMean ); - } else { - - for(;;) { - double u1 = pGenerator -> asDouble(); - double u2 = pGenerator -> asDouble(); - double v1 = 2 * u1 - 1; - double v2 = 2 * u2 - 1; - double w = (v1 * v1) + (v2 * v2); - -// -// We actually generate two IID normal distribution variables. -// We cache the one & return the other. -// - if (w <= 1) { - double y = sqrt( (-2 * log(w)) / w); - double x1 = v1 * y; - double x2 = v2 * y; - - haveCachedNormal = 1; - cachedNormal = x2; - return(x1 * pStdDev + pMean); - } - } - } -} - diff --git a/src/gnu/Normal.h b/src/gnu/Normal.h deleted file mode 100755 index f96ed6dc81..0000000000 --- a/src/gnu/Normal.h +++ /dev/null @@ -1,66 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Normal_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Normal_h - -#include - -class Normal: public Random { - char haveCachedNormal; - double cachedNormal; - -protected: - double pMean; - double pVariance; - double pStdDev; - -public: - Normal(double xmean, double xvariance, RNG *gen); - double mean(); - double mean(double x); - double variance(); - double variance(double x); - virtual double operator()(); -}; - - -inline Normal::Normal(double xmean, double xvariance, RNG *gen) -: Random(gen) { - pMean = xmean; - pVariance = xvariance; - pStdDev = sqrt(pVariance); - haveCachedNormal = 0; -} - -inline double Normal::mean() { return pMean; }; -inline double Normal::mean(double x) { - double t=pMean; pMean = x; - return t; -} - -inline double Normal::variance() { return pVariance; } -inline double Normal::variance(double x) { - double t=pVariance; pVariance = x; - pStdDev = sqrt(pVariance); - return t; -}; - -#endif diff --git a/src/gnu/Poisson.cpp b/src/gnu/Poisson.cpp deleted file mode 100755 index 77a3b96cdf..0000000000 --- a/src/gnu/Poisson.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include <../../nrnconf.h> - -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double Poisson::operator()() -{ - double bound = exp(-1.0 * pMean); - int count = 0; - - for (double product = 1.0; - product >= bound; - product *= pGenerator -> asDouble()) { - count++; - } - return(count - 1); -} diff --git a/src/gnu/Poisson.h b/src/gnu/Poisson.h deleted file mode 100755 index 4bc7054353..0000000000 --- a/src/gnu/Poisson.h +++ /dev/null @@ -1,51 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Poisson_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Poisson_h - -#include - -class Poisson: public Random { -protected: - double pMean; -public: - Poisson(double mean, RNG *gen); - - double mean(); - double mean(double x); - - virtual double operator()(); -}; - - -inline Poisson::Poisson(double mean, RNG *gen) -: Random(gen) { - pMean = mean; -} - -inline double Poisson::mean() { return pMean; } -inline double Poisson::mean(double x) { - double t = pMean; - pMean = x; - return t; -} - -#endif diff --git a/src/gnu/RNG_random123.cpp b/src/gnu/RNG_random123.cpp new file mode 100755 index 0000000000..f4b5711faf --- /dev/null +++ b/src/gnu/RNG_random123.cpp @@ -0,0 +1,12 @@ +#include + +static char initialized = 0; + +RNG_random123::RNG_random123() { + if (!initialized) { + initialized = 1; + c = {{}}; + k = {{}}; + longmurng = std::make_unique>(c.incr(), k); + } +}; diff --git a/src/gnu/RNG_random123.h b/src/gnu/RNG_random123.h new file mode 100755 index 0000000000..74acb05458 --- /dev/null +++ b/src/gnu/RNG_random123.h @@ -0,0 +1,21 @@ +#pragma once + +#include + +#include +#include + +// +// Base class for Random Number Generators using Random123. +// +class RNG_random123 { + r123::Philox4x32::ctr_type c; + r123::Philox4x32::key_type k; + std::shared_ptr> longmurng; + + public: + RNG_random123(); + std::shared_ptr> get_random_generator() { + return longmurng; + } +}; diff --git a/src/gnu/Random.h b/src/gnu/Random.h index a4b9b5c5ec..f86e9f398d 100755 --- a/src/gnu/Random.h +++ b/src/gnu/Random.h @@ -1,20 +1,20 @@ // This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ +/* + Copyright (C) 1988 Free Software Foundation + written by Dirk Grunwald (grunwald@cs.uiuc.edu) + + This file is part of the GNU C++ Library. This library is free + software; you can redistribute it and/or modify it under the terms of + the GNU Library General Public License as published by the Free + Software Foundation; either version 2 of the License, or (at your + option) any later version. This library is distributed in the hope + that it will be useful, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the GNU Library General Public License for more details. + You should have received a copy of the GNU Library General Public + License along with this library; if not, write to the Free Software + Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + */ #ifndef _Random_h #define _Random_h 1 #ifdef __GNUG__ @@ -24,21 +24,23 @@ Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. #include #if MAC - #define Random gnu_Random +#define Random gnu_Random #endif #include +#include class Random { -protected: - RNG *pGenerator; -public: - Random(RNG *generator); + protected: + RNG* pGenerator; + + public: + Random(RNG* generator); virtual ~Random() {} virtual double operator()() = 0; - RNG *generator(); - void generator(RNG *p); + RNG* generator(); + void generator(RNG* p); }; @@ -57,4 +59,183 @@ inline void Random::generator(RNG *p) pGenerator = p; } +class Random_random123 { + protected: + RNG_random123* pGenerator; + + public: + Random_random123(RNG_random123* generator); + virtual ~Random_random123() {} + virtual double operator()() = 0; + + RNG_random123* generator(); + void generator(RNG_random123* p); +}; + + +inline Random_random123::Random_random123(RNG_random123* gen) { + pGenerator = gen; +} + +inline RNG_random123* Random_random123::generator() { + return (pGenerator); +} + +inline void Random_random123::generator(RNG_random123* p) { + pGenerator = p; +} + +#include +class Binomial_random123: public Random_random123 { + protected: + std::binomial_distribution<> d; + + public: + inline Binomial_random123(int t, double p, RNG_random123* gen) + : Random_random123(gen) + , d(t, p) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class DiscreteUniform_random123: public Random_random123 { + protected: + std::uniform_int_distribution<> d; + + public: + inline DiscreteUniform_random123(double low, double high, RNG_random123* gen) + : Random_random123(gen) + , d(low, high) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class Erlang_random123: public Random_random123 { + protected: + double pMean; + double pVariance; + int a; + double b; + inline void setState() { + a = int((pMean * pMean) / pVariance + 0.5); + a = (a > 0) ? a : 1; + b = pMean / a; + } + + std::gamma_distribution<> d; + + public: + inline Erlang_random123(double mean, double variance, RNG_random123* gen) + : Random_random123(gen) + , pMean(mean) + , pVariance(variance) { + setState(); + d = std::gamma_distribution<>(a, b); + } + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class Geometric_random123: public Random_random123 { + protected: + std::geometric_distribution<> d; + + public: + Geometric_random123(double mean, RNG_random123* gen) + : Random_random123(gen) + , d(mean) {} + + double operator()() override { + return (1. + d(*(generator()->get_random_generator()))); + } +}; + +class Normal_random123: public Random_random123 { + protected: + std::normal_distribution<> d; + + public: + inline Normal_random123(double mean, double variance, RNG_random123* gen) + : Random_random123(gen) + , d(mean, std::sqrt(variance)) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class LogNormal_random123: public Normal_random123 { + public: + inline LogNormal_random123(double mean, double variance, RNG_random123* gen) + : Normal_random123(std::log(mean * mean / std::sqrt(variance + (mean * mean))), + log((variance + (mean * mean)) / (mean * mean)), + gen) { + } + + double operator()() override { + return (std::exp(this->Normal_random123::operator()())); + } +}; + +class NegativeExpntl_random123: public Random_random123 { + protected: + std::exponential_distribution<> d; + + public: + inline NegativeExpntl_random123(double xmean, RNG_random123* gen) + : Random_random123(gen) + , d(1 / xmean) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class Poisson_random123: public Random_random123 { + protected: + std::poisson_distribution<> d; + + public: + inline Poisson_random123(double mean, RNG_random123* gen) + : Random_random123(gen) + , d(mean) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class Uniform_random123: public Random_random123 { + protected: + std::uniform_real_distribution<> d; + + public: + inline Uniform_random123(double low, double high, RNG_random123* gen) + : Random_random123(gen) + , d(low, high) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; + +class Weibull_random123: public Random_random123 { + protected: + std::weibull_distribution<> d; + + public: + inline Weibull_random123(double alpha, double beta, RNG_random123* gen) + : Random_random123(gen) + , d(alpha, beta) {} + + inline double operator()() override { + return (d(*(generator()->get_random_generator()))); + } +}; #endif diff --git a/src/gnu/Uniform.cpp b/src/gnu/Uniform.cpp deleted file mode 100755 index ee76bce881..0000000000 --- a/src/gnu/Uniform.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -double Uniform::operator()() -{ - return( pLow + delta * pGenerator -> asDouble() ); -} diff --git a/src/gnu/Uniform.h b/src/gnu/Uniform.h deleted file mode 100755 index c9f3ec63c9..0000000000 --- a/src/gnu/Uniform.h +++ /dev/null @@ -1,71 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Uniform_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Uniform_h 1 - -#include - -// -// The interval [lo..hi] -// - -class Uniform: public Random { - double pLow; - double pHigh; - double delta; -public: - Uniform(double low, double high, RNG *gen); - - double low(); - double low(double x); - double high(); - double high(double x); - - virtual double operator()(); -}; - - -inline Uniform::Uniform(double low, double high, RNG *gen) : Random(gen) -{ - pLow = (low < high) ? low : high; - pHigh = (low < high) ? high : low; - delta = pHigh - pLow; -} - -inline double Uniform::low() { return pLow; } - -inline double Uniform::low(double x) { - double tmp = pLow; - pLow = x; - delta = pHigh - pLow; - return tmp; -} - -inline double Uniform::high() { return pHigh; } - -inline double Uniform::high(double x) { - double tmp = pHigh; - pHigh = x; - delta = pHigh - pLow; - return tmp; -} - -#endif diff --git a/src/gnu/Weibull.cpp b/src/gnu/Weibull.cpp deleted file mode 100755 index 8d1f6b4f74..0000000000 --- a/src/gnu/Weibull.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include <../../nrnconf.h> -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include -#include - -// -// See Simulation, Modelling & Analysis by Law & Kelton, pp259 -// -// This is the ``polar'' method. -// - -double Weibull::operator()() -{ - return( pow(pBeta * ( - log(1 - pGenerator -> asDouble()) ), pInvAlpha) ); -} diff --git a/src/gnu/Weibull.h b/src/gnu/Weibull.h deleted file mode 100755 index 1cc8fd2b13..0000000000 --- a/src/gnu/Weibull.h +++ /dev/null @@ -1,74 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library is distributed in the hope -that it will be useful, but WITHOUT ANY WARRANTY; without even the -implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the GNU Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _Weibull_h -#ifdef __GNUG__ -//#pragma interface -#endif -#define _Weibull_h - -#include - -class Weibull: public Random { -protected: - double pAlpha; - double pInvAlpha; - double pBeta; - - void setState(); - -public: - Weibull(double alpha, double beta, RNG *gen); - - double alpha(); - double alpha(double x); - - double beta(); - double beta(double x); - - virtual double operator()(); -}; - - -inline void Weibull::setState() { - pInvAlpha = 1.0 / pAlpha; -} - -inline Weibull::Weibull(double alpha, double beta, RNG *gen) : Random(gen) -{ - pAlpha = alpha; - pBeta = beta; - setState(); -} - -inline double Weibull::alpha() { return pAlpha; } - -inline double Weibull::alpha(double x) { - double tmp = pAlpha; - pAlpha = x; - setState(); - return tmp; -} - -inline double Weibull::beta() { return pBeta; }; -inline double Weibull::beta(double x) { - double tmp = pBeta; - pBeta = x; - return tmp; -}; - -#endif diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index 27ac546820..3adc50c6c0 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -2,8 +2,8 @@ // definition of random number generation from the g++ library -#include -#include +#include +#include #include "random1.h" #include @@ -17,21 +17,10 @@ #include #include -#include -#include +#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include #include -#include #if HAVE_IV #include "ivoc.h" @@ -202,12 +191,10 @@ void RandomPlay::update(Observable*) { list_remove(); } -Rand::Rand(unsigned long seed, int size, Object* obj) { +Rand::Rand() { // printf("Rand\n"); - gen = new ACG(seed, size); - rand = new Normal(0., 1., gen); - type_ = 0; - obj_ = obj; + gen = new RNG_random123(); + rand = new Normal_random123(0., 1., gen); } Rand::~Rand() { @@ -232,7 +219,7 @@ static void* r_cons(Object* obj) { if (ifarg(2)) size = int(chkarg(2, 7, 98)); - Rand* r = new Rand(seed, size, obj); + Rand* r = new Rand(); return (void*) r; } @@ -242,56 +229,6 @@ static void r_destruct(void* r) { delete (Rand*) r; } -// Use a variant of the Linear Congruential Generator (algorithm M) -// described in Knuth, Art of Computer Programming, Vol. III in -// combination with a Fibonacci Additive Congruential Generator. -// This is a "very high quality" random number generator, -// Default size is 55, giving a size of 1244 bytes to the structure -// minimum size is 7 (total 100 bytes), maximum size is 98 (total 2440 bytes) -// syntax: -// r.ACG([seed],[size]) - -static double r_ACG(void* r) { - Rand* x = (Rand*) r; - - unsigned long seed = 0; - int size = 55; - - if (ifarg(1)) - seed = long(*getarg(1)); - if (ifarg(2)) - size = int(chkarg(2, 7, 98)); - - x->rand->generator(new ACG(seed, size)); - x->type_ = 0; - delete x->gen; - x->gen = x->rand->generator(); - return 1.; -} - -// Use a Multiplicative Linear Congruential Generator. Not as high -// quality as the ACG, but uses only 8 bytes -// syntax: -// r.MLCG([seed1],[seed2]) - -static double r_MLCG(void* r) { - Rand* x = (Rand*) r; - - unsigned long seed1 = 0; - unsigned long seed2 = 0; - - if (ifarg(1)) - seed1 = long(*getarg(1)); - if (ifarg(2)) - seed2 = long(*getarg(2)); - - x->rand->generator(new MLCG(seed1, seed2)); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 1; - return 1.; -} - static double r_MCellRan4(void* r) { Rand* x = (Rand*) r; @@ -425,6 +362,11 @@ static double r_repick(void* r) { return (*(x->rand))(); } +static double r_repick_random123(void* r) { + Rand* x = (Rand*) r; + return (*(x->rand))(); +} + double nrn_random_pick(Rand* r) { if (r) { return (*(r->rand))(); @@ -449,103 +391,90 @@ Rand* nrn_random_arg(int i) { // uniform random variable over the open interval [low...high) // syntax: // r.uniform(low,high) - static double r_uniform(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); double a2 = *getarg(2); delete x->rand; - x->rand = new Uniform(a1, a2, x->gen); + x->rand = new Uniform_random123(a1, a2, x->gen); return (*(x->rand))(); } // uniform random variable over the closed interval [low...high] // syntax: // r.discunif(low,high) - static double r_discunif(void* r) { Rand* x = (Rand*) r; long a1 = long(*getarg(1)); long a2 = long(*getarg(2)); delete x->rand; - x->rand = new DiscreteUniform(a1, a2, x->gen); + x->rand = new DiscreteUniform_random123(a1, a2, x->gen); return (*(x->rand))(); } - // normal (gaussian) distribution // syntax: // r.normal(mean,variance) - static double r_normal(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); double a2 = *getarg(2); delete x->rand; - x->rand = new Normal(a1, a2, x->gen); + x->rand = new Normal_random123(a1, a2, x->gen); return (*(x->rand))(); } - // logarithmic normal distribution // syntax: // r.lognormal(mean) - static double r_lognormal(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); double a2 = *getarg(2); delete x->rand; - x->rand = new LogNormal(a1, a2, x->gen); + x->rand = new LogNormal_random123(a1, a2, x->gen); return (*(x->rand))(); } - // poisson distribution // syntax: // r.poisson(mean) - static double r_poisson(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); delete x->rand; - x->rand = new Poisson(a1, x->gen); + x->rand = new Poisson_random123(a1, x->gen); return (*(x->rand))(); } - // binomial distribution, which models successfully drawing items from a pool // n is the number items in the pool and p is the probablity of each item // being successfully drawn (n>0, 0<=p<=1) // syntax: // r.binomial(n,p) - static double r_binomial(void* r) { Rand* x = (Rand*) r; int a1 = int(chkarg(1, 0, 1e99)); double a2 = chkarg(2, 0, 1); delete x->rand; - x->rand = new Binomial(a1, a2, x->gen); + x->rand = new Binomial_random123(a1, a2, x->gen); return (*(x->rand))(); } - // discrete geometric distribution // Given 0<=mean<=1, returns the number of uniform random samples // that were drawn before the sample was larger than mean (always // greater than 0. // syntax: // r.geometric(mean) - static double r_geometric(void* r) { Rand* x = (Rand*) r; double a1 = chkarg(1, 0, 1); delete x->rand; - x->rand = new Geometric(a1, x->gen); + x->rand = new Geometric_random123(a1, x->gen); return (*(x->rand))(); } - // hypergeometric distribution // syntax: // r.hypergeo(mean,variance) @@ -563,26 +492,23 @@ static double r_hypergeo(void* r) { // negative exponential distribution // syntax: // r.negexp(mean) - static double r_negexp(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); delete x->rand; - x->rand = new NegativeExpntl(a1, x->gen); + x->rand = new NegativeExpntl_random123(a1, x->gen); return (*(x->rand))(); } - // Erlang distribution // syntax: // r.erlang(mean,variance) - static double r_erlang(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); double a2 = *getarg(2); delete x->rand; - x->rand = new Erlang(a1, a2, x->gen); + x->rand = new Erlang_random123(a1, a2, x->gen); return (*(x->rand))(); } @@ -590,13 +516,12 @@ static double r_erlang(void* r) { // Weibull distribution // syntax: // r.weibull(alpha,beta) - static double r_weibull(void* r) { Rand* x = (Rand*) r; double a1 = *getarg(1); double a2 = *getarg(2); delete x->rand; - x->rand = new Weibull(a1, a2, x->gen); + x->rand = new Weibull_random123(a1, a2, x->gen); return (*(x->rand))(); } @@ -612,14 +537,13 @@ extern "C" void nrn_random_play() { } -static Member_func r_members[] = {{"ACG", r_ACG}, - {"MLCG", r_MLCG}, - {"Isaac64", r_Isaac64}, +static Member_func r_members[] = {{"Isaac64", r_Isaac64}, {"MCellRan4", r_MCellRan4}, {"Random123", r_nrnran123}, {"Random123_globalindex", r_ran123_globalindex}, {"seq", r_sequence}, {"repick", r_repick}, + {"repick_random123", r_repick_random123}, {"uniform", r_uniform}, {"discunif", r_discunif}, {"normal", r_normal}, @@ -660,4 +584,4 @@ int nrn_random123_setseq(void* r, uint32_t seq, char which) { } void nrn_set_random_sequence(void* r, int seq) { nrn_set_random_sequence(static_cast(r), static_cast(seq)); -} \ No newline at end of file +} diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index ab10014186..ecddf85214 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -123,7 +123,6 @@ static double dmaxint_ = 9007199254740992; // definition of random numer generator #include "random1.h" -#include #if HAVE_IV #include "utility.h" diff --git a/src/ivoc/random1.h b/src/ivoc/random1.h index 8fedd8f5b0..6952a414c6 100644 --- a/src/ivoc/random1.h +++ b/src/ivoc/random1.h @@ -1,5 +1,4 @@ -#ifndef random1_h -#define random1_h +#pragma once #include "RNG.h" #include "Random.h" @@ -8,13 +7,11 @@ struct Object; class Rand { public: - Rand(unsigned long seed = 0, int size = 55, Object* obj = NULL); + Rand(); ~Rand(); - RNG* gen; - Random* rand; + RNG_random123* gen; + Random_random123* rand; int type_; // can do special things with some kinds of RNG // double* looks like random variable that gets changed on every fadvance Object* obj_; }; - -#endif diff --git a/src/nrniv/singlech.cpp b/src/nrniv/singlech.cpp index 3746d13fc5..e7efbd5811 100644 --- a/src/nrniv/singlech.cpp +++ b/src/nrniv/singlech.cpp @@ -8,7 +8,6 @@ #include "singlech.h" #include "membfunc.h" #include "random1.h" -#include "NegExp.h" extern "C" { diff --git a/test_random.py b/test_random.py new file mode 100644 index 0000000000..8752fd649c --- /dev/null +++ b/test_random.py @@ -0,0 +1,94 @@ +from neuron import h +import matplotlib.pyplot as plt + + +def create_plot(rand, type): + values = [] + for i in range(100000): + values.append(rand.repick()) + # plt.plot(values) + plt.hist(values, bins=1000) + # plt.show() + plt.savefig("hist_{}.pdf".format(type), format="pdf") + plt.close() + + +rand = h.Random() +print("Test default rand") +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "random") + +print("Setup rand for Binomial distribution") +print(rand.binomial(20, 0.5)) +# print(rand.repick()) +# print(rand.binomial(20, 0.5)) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "binomial") + +print("Setup rand for Normal distribution") +print(rand.normal(5, 2)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "normal") + +print("Setup rand for LogNormal distribution") +print(rand.lognormal(5, 2)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "lognormal") + +print("Setup rand for Poisson distribution") +print(rand.poisson(5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "poisson") + +print("Setup rand for Uniform distribution") +print(rand.uniform(1, 5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "uniform") + +print("Setup rand for DiscUnif distribution") +print(rand.discunif(1, 5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "discunif") + +print("Setup rand for Weibull distribution") +print(rand.weibull(1, 5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "weibull") + +print("Setup rand for Erlang distribution") +print(rand.erlang(1, 5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "erlang") + +print("Setup rand for Erlang distribution") +print(rand.negexp(1.5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "negexp") + +print("Setup rand for Geom distribution") +print(rand.geometric(0.5)) +# print(rand.repick()) +# for i in range(10): +# print(rand.repick()) +create_plot(rand, "geom") + +quit() diff --git a/test_random123.py b/test_random123.py new file mode 100644 index 0000000000..010a257e63 --- /dev/null +++ b/test_random123.py @@ -0,0 +1,94 @@ +from neuron import h +import matplotlib.pyplot as plt + + +def create_plot(rand, type): + values = [] + for i in range(100000): + values.append(rand.repick_random123()) + # plt.plot(values) + plt.hist(values, bins=1000) + # plt.show() + plt.savefig("hist_{}.pdf".format(type), format="pdf") + plt.close() + + +rand = h.Random() +print("Test default rand") +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "random") + +print("Setup rand for Binomial distribution") +print(rand.binomial_random123(20, 0.5)) +# print(rand.repick_random123()) +# print(rand.binomial_random123(20, 0.5)) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "binomial") + +print("Setup rand for Normal distribution") +print(rand.normal_random123(5, 2)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "normal") + +print("Setup rand for LogNormal distribution") +print(rand.lognormal_random123(5, 2)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "lognormal") + +print("Setup rand for Poisson distribution") +print(rand.poisson_random123(5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "poisson") + +print("Setup rand for Uniform distribution") +print(rand.uniform_random123(1, 5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "uniform") + +print("Setup rand for DiscUnif distribution") +print(rand.discunif_random123(1, 5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "discunif") + +print("Setup rand for Weibull distribution") +print(rand.weibull_random123(1, 5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "weibull") + +print("Setup rand for Erlang distribution") +print(rand.erlang_random123(1, 5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "erlang") + +print("Setup rand for NegExp distribution") +print(rand.negexp_random123(1.5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "negexp") + +print("Setup rand for Geom distribution") +print(rand.geometric_random123(0.5)) +# print(rand.repick_random123()) +# for i in range(10): +# print(rand.repick_random123()) +create_plot(rand, "geom") + +quit()