Commit 699c7722 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

[FIX] (libtsxx): fix random number sequence

Multiple calls to ts::rnd::dugauss() resulted in identical sequences of random
numbers.

Issue regarding random number generator in libtsxx
==================================================

Location: src/libs/libtsxx/random.cc

The function ts::rnd::dugauss() allocates a GSL random number generator and
immediately seeds it. The seeding operation is required, since the default
seed is a constant, which would result in identical sequences. Function
ts::rnd::dugauss() uses the current system time as returned by function time()
as a seed value. Unfortunately, if ts::rnd::dugauss() is called several times
during the same second, this also will result in identical sequences.

Currently ts::rnd::dugauss() is only used by src/ts/wf/noisymize.cc and by the
filter GaussianNoise in libtsxx.

The appropriate solution appears to be:
- encapsulate the GSL randon number generator in a class object, just like is
  provided in tfxx::numeric::RNGgaussian
- use a static variable of this class type in ts::rnd::dugauss()
  this variable will be initialized only once upon the first call to the
  function
- the destructor is guaranteed to be called after termination of the program,
  such that the destructor is able to properly call gsl_rng_free
parent adc184bb
......@@ -4,11 +4,11 @@
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 27/06/2006
* \date 19/03/2016
*
* create a random series (implementation)
*
* Copyright (c) 2006 by Thomas Forbriger (BFO Schiltach)
* Copyright (c) 2006, 2016 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
......@@ -29,11 +29,12 @@
*
* REVISIONS and CHANGES
* - 27/06/2006 V1.0 Thomas Forbriger
* - 19/03/2016 V1.1 seed sequence only once
*
* ============================================================================
*/
#define TS_RANDOM_CC_VERSION \
"TS_RANDOM_CC V1.0"
"TS_RANDOM_CC V1.1"
#include <tsxx/random.h>
#include <gsl/gsl_rng.h>
......@@ -44,20 +45,39 @@ namespace ts {
namespace rnd {
namespace helper {
//! internal class used for static variable in dugauss
class GSLrng {
public:
explicit GSLrng(const gsl_rng_type* T)
{
gsl_rng_env_setup();
MR=gsl_rng_alloc(T);
this->seed();
}
~GSLrng() { gsl_rng_free(MR); }
double ugaussian() const { return gsl_ran_ugaussian(MR); }
void seed(const unsigned long int& seedvalue) const
{ gsl_rng_set(MR, seedvalue); }
void seed() const { this->seed(time(0)); }
private:
gsl_rng* MR;
}; // class GSLrng
} // namespace helper
/* ---------------------------------------------------------------------- */
Tdseries dugauss(const int& n)
{
static helper::GSLrng rng(gsl_rng_default);
Tdseries retval(n);
gsl_rng_env_setup();
const gsl_rng_type* T=gsl_rng_default;
gsl_rng* R=gsl_rng_alloc(T);
gsl_rng_set(R, time(0));
for (int i=retval.f(); i<=retval.l(); ++i)
{ retval(i)=gsl_ran_ugaussian(R); }
gsl_rng_free(R);
{ retval(i)=rng.ugaussian(); }
return (retval);
} // Tdseries dugauss
} // namespace rnd
} // namespace ts
......
......@@ -4,7 +4,7 @@
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
* \date 20/12/2003
* \date 19/03/2016
*
* test time series modules
*
......@@ -33,11 +33,12 @@
* - 21/11/2006 V1.2 test drop containers
* - 28/01/2012 V1.3 test OffsetVariableTaper
* - 20/05/2015 V1.4 print library usage texts
* - 19/03/2016 V1.5 test subsequent calls to function dugauss
*
* ============================================================================
*/
#define TSTEST_VERSION \
"TSTEST V1.4 test time series modules"
"TSTEST V1.5 test time series modules"
#include <tfxx/commandline.h>
#include <iostream>
......@@ -182,7 +183,9 @@ int main(int iargc, char* argv[])
if (opt.randomnoise)
{
cout << "test random noise generation" << endl;
ts::rnd::Tdseries s=ts::rnd::dugauss(30);
ts::rnd::Tdseries s=ts::rnd::dugauss(10);
DUMP( s );
s=ts::rnd::dugauss(10);
DUMP( s );
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment