portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2024-11-10 00:45:36 +00:00
Code Releases Activity

internal random generator entirely removed, use C++11 generators

Browse Source
This commit is contained in:
Antonin Portelli 2016-04-06 20:04:23 +01:00
parent 984947c677
commit 3611d12fa9
20 changed files with 71 additions and 1044 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
examples/exFit.cpp
View File

@ -20,11 +20,13 @@ const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void)
{
// generate fake data
XYStatData data;
RandGen rg;
double xBuf[2];
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
XYStatData data;
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
double xBuf[2];
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
cout << "-- generating fake data..." << endl;
data.addXDim(nPoint1);
@ -33,13 +35,13 @@ int main(void)
for (Index i1 = 0; i1 < nPoint1; ++i1)
{
xBuf[0] = i1*dx1;
data.x(i1, 0) = rg.gaussian(xBuf[0], xErr);
data.x(i1, 0) = xErr*dis(gen) + xBuf[0];
for (Index i2 = 0; i2 < nPoint2; ++i2)
{
xBuf[1] = i2*dx2;
data.x(i2, 1) = xBuf[1];
data.y(data.dataIndex(i1, i2), 0) = rg.gaussian(f(xBuf, exactPar),
yErr);
data.y(data.dataIndex(i1, i2), 0) = yErr*dis(gen)
+ f(xBuf, exactPar);
}
}
data.setXError(0, DVec::Constant(data.getXSize(0), xErr));

18
examples/exFitSample.cpp
View File

@ -20,11 +20,13 @@ const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void)
{
// generate fake data
XYSampleData data(nSample);
RandGen rg;
double xBuf[2];
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
XYSampleData data(nSample);
double xBuf[2];
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
cout << "-- generating fake data..." << endl;
data.addXDim(nPoint1);
@ -35,13 +37,13 @@ int main(void)
for (Index i1 = 0; i1 < nPoint1; ++i1)
{
xBuf[0] = i1*dx1;
data.x(i1, 0)[s] = rg.gaussian(xBuf[0], xErr);
data.x(i1, 0)[s] = xErr*dis(gen) + xBuf[0];
for (Index i2 = 0; i2 < nPoint2; ++i2)
{
xBuf[1] = i2*dx2;
data.x(i2, 1)[s] = xBuf[1];
data.y(data.dataIndex(i1, i2), 0)[s] =
rg.gaussian(f(xBuf, exactPar), yErr);
data.y(data.dataIndex(i1, i2), 0)[s] = yErr*dis(gen)
+ f(xBuf, exactPar);
}
}
}

40
examples/exRand.cpp
View File

@ -2,52 +2,26 @@
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/Plot.hpp>
#include <LatAnalyze/RandGen.hpp>
using namespace std;
using namespace Latan;
constexpr int seqLength = 25;
constexpr int saveStep = 9;
constexpr Index nDraw = 20000;
const string stateFileName = "exRand.seed";
int main(void)
{
RandGenState state;
RandGen gen[2];
AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read);
DVec gauss(nDraw);
Plot p;
Histogram h;
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
DVec gauss(nDraw);
Plot p;
Histogram h;
cout << "- GENERATOR STATE I/O TESTS" << endl;
cout << "-- generating a " << seqLength << " steps random sequence..."
<< endl;
for (int i = 0; i < seqLength; ++i)
{
if (i == saveStep)
{
state = gen[0].getState();
stateFile.save(state, "exRand");
cout << "generator state after step " << saveStep - 1
<< " saved in '" << stateFileName << "'" << endl;
}
cout << "step " << i << "\t: " << gen[0].uniform() <<endl;
}
cout << "-- setting up another generator from '" << stateFileName << "'..."
<< endl;
gen[1].setState(stateFile.read<RandGenState>("exRand"));
cout << "-- generating a " << seqLength << " steps random sequence..."
<< endl;
for (int i = 0; i < seqLength; ++i)
{
cout << "step " << i << "\t: " << gen[1].uniform() <<endl;
}
cout << "-- generating " << nDraw << " Gaussian random numbers..." << endl;
FOR_VEC(gauss, i)
{
gauss(i) = gen[0].gaussian();
gauss(i) = dis(gen);
}
h.setFromData(gauss, -5., 5., 40);
h.normalize();

14
lib/AsciiFile.cpp
View File

@ -105,20 +105,6 @@ void AsciiFile::save(const DMatSample &ms, const std::string &name)
fileStream_ << "#L latan_end rs_sample " << endl;
}
void AsciiFile::save(const RandGenState &state, const std::string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin rg_state " << name << endl;
fileStream_ << state << endl;
fileStream_ << "#L latan_end rg_state " << endl;
}
// read first name ////////////////////////////////////////////////////////////
string AsciiFile::getFirstName(void)
{

5
lib/AsciiFile.hpp
View File

@ -25,7 +25,6 @@
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/ParserState.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
BEGIN_LATAN_NAMESPACE
@ -48,12 +47,11 @@ public:
bool isFirst;
std::string first;
// parsing buffers
RandGenState stateBuf;
int intBuf;
DSample dSampleBuf;
DMatSample dMatSampleBuf;
std::queue<DMat> dMatQueue;
std::queue<double> doubleQueue;
std::queue<int> intQueue;
private:
// allocation/deallocation functions defined in IoAsciiLexer.lpp
virtual void initScanner(void);
@ -69,7 +67,6 @@ public:
virtual void save(const DMat &m, const std::string &name);
virtual void save(const DSample &ds, const std::string &name);
virtual void save(const DMatSample &ms, const std::string &name);
virtual void save(const RandGenState &state, const std::string &name);
// read first name
virtual std::string getFirstName(void);
// tests

2
lib/AsciiLexer.lpp
View File

@ -76,7 +76,7 @@ BLANK [ \t]
%%
{LMARK} {BEGIN(MARK);}
{INT} {yylval->val_int = strTo<int>(yytext); RETTOK(INT);}
{INT} {yylval->val_int = strTo<long int>(yytext); RETTOK(INT);}
{FLOAT} {yylval->val_double = strTo<double>(yytext); RETTOK(FLOAT);}
({ALPHA}|{DIGIT})+ {strcpy(yylval->val_str,yytext); RETTOK(ID);}
<MARK>latan_begin {BEGIN(TYPE); RETTOK(OPEN);}

44
lib/AsciiParser.ypp
View File

@ -22,11 +22,6 @@
#include <LatAnalyze/AsciiFile.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <iostream>
#include <sstream>
#include <utility>
#include <cstring>
#if (defined __INTEL_COMPILER)
#pragma warning disable 1682 2259
@ -47,7 +42,7 @@
%}
%pure-parser
%name-prefix = "_Ascii_"
%name-prefix "_Ascii_"
%locations
%defines
%error-verbose
@ -57,10 +52,10 @@
%union
{
int val_int;
double val_double;
char val_char;
char val_str[256];
long int val_int;
double val_double;
char val_char;
char val_str[256];
}
%token <val_char> ERR
@ -69,7 +64,7 @@
%token <val_str> ID
%token OPEN CLOSE MAT SAMPLE RG_STATE
%type <val_str> mat matsample dsample rg_state
%type <val_str> mat matsample dsample
%{
int _Ascii_lex(YYSTYPE* lvalp, YYLTYPE* llocp, void* scanner);
@ -111,11 +106,6 @@ data:
TEST_FIRST($1);
(*state->data)[$1].reset(new DMatSample(state->dMatSampleBuf));
}
| rg_state
{
TEST_FIRST($1);
(*state->data)[$1].reset(new RandGenState(state->stateBuf));
}
;
mat:
@ -187,23 +177,6 @@ matsample:
}
;
rg_state:
OPEN RG_STATE ID ints CLOSE RG_STATE
{
if (state->intQueue.size() != RLXG_STATE_SIZE)
{
LATAN_ERROR(Size, "random generator state '" + *state->streamName
+ ":" + $3 + "' has a wrong size");
}
for (Index i = 0; i < RLXG_STATE_SIZE; ++i)
{
state->stateBuf[i] = state->intQueue.front();
state->intQueue.pop();
}
strcpy($$, $3);
}
;
mats:
mats mat
| mat
@ -215,8 +188,3 @@ floats:
| FLOAT {state->doubleQueue.push($1);}
| INT {state->doubleQueue.push(static_cast<double>($1));}
;
ints:
ints INT {state->intQueue.push($2);}
| INT {state->intQueue.push($1);}
;

35
lib/Dataset.hpp
View File

@ -23,9 +23,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/StatArray.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
#include <vector>
BEGIN_LATAN_NAMESPACE
@ -35,6 +32,8 @@ BEGIN_LATAN_NAMESPACE
template <typename T>
class Dataset: public StatArray<T>
{
public:
typedef std::random_device::result_type SeedType;
public:
// constructors
Dataset(void) = default;
@ -46,7 +45,8 @@ public:
template <typename FileType>
void load(const std::string &listFileName, const std::string &dataName);
// resampling
Sample<T> bootstrapMean(const Index nSample, RandGen& generator);
Sample<T> bootstrapMean(const Index nSample, const SeedType seed);
Sample<T> bootstrapMean(const Index nSample);
private:
// mean from pointer vector for resampling
void ptVectorMean(T &m, const std::vector<const T *> &v);
@ -82,26 +82,23 @@ void Dataset<T>::load(const std::string &listFileName,
// resampling //////////////////////////////////////////////////////////////////
template <typename T>
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, RandGen& generator)
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, const SeedType seed)
{
typedef typename std::vector<const T *>::size_type size_type;
size_type nData = static_cast<size_type>(this->size());
std::vector<const T *> data(nData);
Sample<T> s(nSample);
std::vector<const T *> data(this->size());
Sample<T> s(nSample);
std::mt19937 gen(seed);
std::uniform_int_distribution<Index> dis(0, this->size() - 1);
for (size_type j = 0; j < nData; ++j)
for (unsigned int j = 0; j < this->size(); ++j)
{
data[j] = &((*this)[static_cast<Index>(j)]);
}
ptVectorMean(s[central], data);
for (Index i = 0; i < nSample; ++i)
{
for (size_type j = 0; j < nData; ++j)
for (unsigned int j = 0; j < this->size(); ++j)
{
Index k= static_cast<Index>(generator.discreteUniform(static_cast<unsigned int>(nData)));
data[j] = &((*this)[k]);
data[j] = &((*this)[dis(gen)]);
}
ptVectorMean(s[i], data);
}
@ -109,6 +106,14 @@ Sample<T> Dataset<T>::bootstrapMean(const Index nSample, RandGen& generator)
return s;
}
template <typename T>
Sample<T> Dataset<T>::bootstrapMean(const Index nSample)
{
std::random_device rd;
return bootstrapMean(nSample, rd());
}
template <typename T>
void Dataset<T>::ptVectorMean(T &m, const std::vector<const T *> &v)
{

4
lib/File.hpp
View File

@ -24,9 +24,6 @@
#include <LatAnalyze/IoObject.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <queue>
#include <unordered_map>
BEGIN_LATAN_NAMESPACE
@ -63,7 +60,6 @@ public:
virtual void save(const DMat &m, const std::string &name) = 0;
virtual void save(const DSample &ds, const std::string &name) = 0;
virtual void save(const DMatSample &ms, const std::string &name) = 0;
virtual void save(const RandGenState &state, const std::string &name) = 0;
// read first name
virtual std::string getFirstName(void) = 0;
// tests

45
lib/Hdf5File.cpp
View File

@ -32,7 +32,6 @@ constexpr unsigned int maxGroupNameSize = 1024u;
const short dMatType = static_cast<short>(IoObject::IoType::dMat);
const short dSampleType = static_cast<short>(IoObject::IoType::dSample);
const short dMatSampleType = static_cast<short>(IoObject::IoType::dMatSample);
const short rgStateType = static_cast<short>(IoObject::IoType::rgState);
/******************************************************************************
* Hdf5File implementation *
@ -131,27 +130,6 @@ void Hdf5File::save(const DMatSample &ms, const string &name)
}
}
void Hdf5File::save(const RandGenState &state, const string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
Group group;
Attribute attr;
DataSet dataset;
hsize_t dim = RLXG_STATE_SIZE;
hsize_t attrDim = 1;
DataSpace dataSpace(1, &dim), attrSpace(1, &attrDim);
group = h5File_->createGroup(name.c_str() + nameOffset(name));
attr = group.createAttribute("type", PredType::NATIVE_SHORT, attrSpace);
attr.write(PredType::NATIVE_SHORT, &rgStateType);
dataset = group.createDataSet("data", PredType::NATIVE_INT, dataSpace);
dataset.write(state.data(), PredType::NATIVE_INT);
}
// read first name ////////////////////////////////////////////////////////////
string Hdf5File::getFirstName(void)
{
@ -288,20 +266,6 @@ void Hdf5File::load(DSample &ds, const DataSet &d)
d.read(ds.data(), PredType::NATIVE_DOUBLE);
}
void Hdf5File::load(RandGenState &state, const DataSet &d)
{
DataSpace dataspace;
hsize_t dim[1];
dataspace = d.getSpace();
dataspace.getSimpleExtentDims(dim);
if (dim[0] != RLXG_STATE_SIZE)
{
LATAN_ERROR(Io, "random generator state has a wrong length");
}
d.read(state.data(), PredType::NATIVE_INT);
}
string Hdf5File::load(const string &name)
{
if ((mode_ & Mode::read)&&(isOpen()))
@ -364,15 +328,6 @@ string Hdf5File::load(const string &name)
}
break;
}
case IoObject::IoType::rgState:
{
RandGenState *pt = new RandGenState;
data_[groupName].reset(pt);
dataset = group.openDataSet("data");
load(*pt, dataset);
break;
}
default:
{
LATAN_ERROR(Io, "unknown data type ("

3
lib/Hdf5File.hpp
View File

@ -24,7 +24,6 @@
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <H5Cpp.h>
BEGIN_LATAN_NAMESPACE
@ -48,7 +47,6 @@ public:
virtual void save(const DMat &m, const std::string &name);
virtual void save(const DSample &ds, const std::string &name);
virtual void save(const DMatSample &ms, const std::string &name);
virtual void save(const RandGenState &state, const std::string &name);
// read first name
virtual std::string getFirstName(void);
// tests
@ -63,7 +61,6 @@ private:
void load(DMat &m, const H5NS::DataSet &d);
void load(DSample &ds, const H5NS::DataSet &d);
void load(DMatSample &s, const H5NS::DataSet &d);
void load(RandGenState &state, const H5NS::DataSet &d);
// check name for forbidden characters
static size_t nameOffset(const std::string &name);
private:

3
lib/IoObject.hpp
View File

@ -36,8 +36,7 @@ public:
noType = 0,
dMat = 1,
dMatSample = 2,
rgState = 3,
dSample = 4
dSample = 3
};
public:
// destructor

2
lib/Makefile.am
View File

@ -40,7 +40,6 @@ libLatAnalyze_la_SOURCES = \
Minimizer.cpp \
Model.cpp \
Plot.cpp \
RandGen.cpp \
RootFinder.cpp \
Solver.cpp \
StatArray.cpp \
@ -75,7 +74,6 @@ libLatAnalyze_la_HEADERS = \
Model.hpp \
ParserState.hpp \
Plot.hpp \
RandGen.hpp \
RootFinder.hpp \
TabFunction.hpp \
Solver.hpp \

5
lib/MathParser.ypp
View File

@ -18,9 +18,6 @@
*/
%{
#include <iostream>
#include <sstream>
#include <cstring>
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
@ -36,7 +33,7 @@
%}
%pure-parser
%name-prefix = "_math_"
%name-prefix "_math_"
%locations
%defines
%error-verbose

690
lib/RandGen.cpp
View File

@ -1,690 +0,0 @@
/*
* RandGen.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/RandGen.hpp>
#include <LatAnalyze/includes.hpp>
#pragma GCC diagnostic ignored "-Wfloat-conversion"
#ifndef RLXD_LEVEL
#define RLXD_LEVEL 1
#endif
using namespace std;
using namespace Latan;
/******************************************************************************
* RandGenState implementation *
******************************************************************************/
// IO type ///////////////////////////////////////////////////////////////
IoObject::IoType RandGenState::getType(void) const
{
return IoType::rgState;
}
/******************************************************************************
* RandGen implementation *
******************************************************************************/
// RanLxd implementation ///////////////////////////////////////////////////////
RandGen::RanLxd::RanLxd(void)
{
// avoid a warning in the SSE case
one_bit = 0.0;
}
/****************************************************************************/
/* Copyright (C) 2005 Martin Luescher (GPL)
* This software is distributed under the terms of the GNU General Public
* License (GPL)
*
* Random number generator "ranlxd". See the notes
*
* "User's guide for ranlxs and ranlxd v3.2" (December 2005)
*
* "Algorithms used in ranlux v3.0" (May 2001)
*
* for a detailed description
*
* The functions are
*
* void ranlxd(double r[],int n)
* Computes the next n double-precision random numbers and
* assigns them to the elements r[0],...,r[n-1] of the array r[]
*
* void rlxd_init(int level,int seed)
* Initialization of the generator
*
* int rlxd_size(void)
* Returns the number of integers required to save the state of
* the generator
*
* void rlxd_get(int state[])
* Extracts the current state of the generator and stores the
* information in the array state[N] where N>=rlxd_size()
*
* void rlxd_reset(int state[])
* Resets the generator to the state defined by the array state[N]
*/
#ifdef HAVE_SSE
#define STEP(pi,pj) \
__asm__ __volatile__ ("movaps %4, %%xmm4 \n\t" \
"movaps %%xmm2, %%xmm3 \n\t" \
"subps %2, %%xmm4 \n\t" \
"movaps %%xmm1, %%xmm5 \n\t" \
"cmpps $0x6, %%xmm4, %%xmm2 \n\t" \
"andps %%xmm2, %%xmm5 \n\t" \
"subps %%xmm3, %%xmm4 \n\t" \
"andps %%xmm0, %%xmm2 \n\t" \
"addps %%xmm4, %%xmm5 \n\t" \
"movaps %%xmm5, %0 \n\t" \
"movaps %5, %%xmm6 \n\t" \
"movaps %%xmm2, %%xmm3 \n\t" \
"subps %3, %%xmm6 \n\t" \
"movaps %%xmm1, %%xmm7 \n\t" \
"cmpps $0x6, %%xmm6, %%xmm2 \n\t" \
"andps %%xmm2, %%xmm7 \n\t" \
"subps %%xmm3, %%xmm6 \n\t" \
"andps %%xmm0, %%xmm2 \n\t" \
"addps %%xmm6, %%xmm7 \n\t" \
"movaps %%xmm7, %1" \
: \
"=m" ((*pi).c1), \
"=m" ((*pi).c2) \
: \
"m" ((*pi).c1), \
"m" ((*pi).c2), \
"m" ((*pj).c1), \
"m" ((*pj).c2) \
: \
"xmm2", "xmm3", "xmm4", "xmm5", "xmm6", "xmm7")
void RandGen::RanLxd::error(int no) const
{
switch(no)
{
case 1:
LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)");
break;
case 2:
LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)");
break;
case 3:
LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)");
break;
case 5:
LATAN_ERROR(Logic, "Unexpected input data");
break;
}
}
void RandGen::RanLxd::update(void)
{
int k,kmax;
rlxd_dble_vec_t *pmin,*pmax,*pi,*pj;
kmax=rlxd_pr;
pmin=&rlxd_x.vec[0];
pmax=pmin+12;
pi=&rlxd_x.vec[ir];
pj=&rlxd_x.vec[jr];
__asm__ __volatile__ ("movaps %0, %%xmm0 \n\t"
"movaps %1, %%xmm1 \n\t"
"movaps %2, %%xmm2"
:
:
"m" (one_bit_sse),
"m" (one_sse),
"m" (carry)
:
"xmm0", "xmm1", "xmm2");
for (k=0;k<kmax;k++)
{
STEP(pi,pj);
pi+=1;
pj+=1;
if (pi==pmax)
pi=pmin;
if (pj==pmax)
pj=pmin;
}
__asm__ __volatile__ ("movaps %%xmm2, %0"
:
"=m" (carry));
ir+=prm;
jr+=prm;
if (ir>=12)
ir-=12;
if (jr>=12)
jr-=12;
is=8*ir;
is_old=is;
}
void RandGen::RanLxd::define_constants(void)
{
int k;
float b;
one_sse.c1=1.0f;
one_sse.c2=1.0f;
one_sse.c3=1.0f;
one_sse.c4=1.0f;
b=(float)(ldexp(1.0,-24));
one_bit_sse.c1=b;
one_bit_sse.c2=b;
one_bit_sse.c3=b;
one_bit_sse.c4=b;
for (k=0;k<96;k++)
{
next[k]=(k+1)%96;
if ((k%4)==3)
next[k]=(k+5)%96;
}
}
void RandGen::RanLxd::rlxd_init(int level,int seed)
{
int i,k,l;
int ibit,jbit,xbit[31];
int ix,iy;
define_constants();
if (level==1)
rlxd_pr=202;
else if (level==2)
rlxd_pr=397;
else
error(1);
i=seed;
for (k=0;k<31;k++)
{
xbit[k]=i%2;
i/=2;
}
if ((seed<=0)||(i!=0))
error(2);
ibit=0;
jbit=18;
for (i=0;i<4;i++)
{
for (k=0;k<24;k++)
{
ix=0;
for (l=0;l<24;l++)
{
iy=xbit[ibit];
ix=2*ix+iy;
xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
ibit=(ibit+1)%31;
jbit=(jbit+1)%31;
}
if ((k%4)!=i)
ix=16777215-ix;
rlxd_x.num[4*k+i]=(float)(ldexp((double)(ix),-24));
}
}
carry.c1=0.0f;
carry.c2=0.0f;
carry.c3=0.0f;
carry.c4=0.0f;
ir=0;
jr=7;
is=91;
is_old=0;
prm=rlxd_pr%12;
init=1;
}
void RandGen::RanLxd::ranlxd(double r[],int n)
{
int k;
if (init==0)
rlxd_init(1,1);
for (k=0;k<n;k++)
{
is=next[is];
if (is==is_old)
update();
r[k]=(double)(rlxd_x.num[is+4])+(double)(one_bit_sse.c1*rlxd_x.num[is]);
}
}
int RandGen::RanLxd::rlxd_size(void) const
{
return(RLXG_STATE_SIZE);
}
void RandGen::RanLxd::rlxd_get(int state[]) const
{
int k;
float base;
if (init==0)
error(3);
base=(float)(ldexp(1.0,24));
state[0]=rlxd_size();
for (k=0;k<96;k++)
state[k+1]=(int)(base*rlxd_x.num[k]);
state[97]=(int)(base*carry.c1);
state[98]=(int)(base*carry.c2);
state[99]=(int)(base*carry.c3);
state[100]=(int)(base*carry.c4);
state[101]=rlxd_pr;
state[102]=ir;
state[103]=jr;
state[104]=is;
}
void RandGen::RanLxd::rlxd_reset(const int state[])
{
int k;
define_constants();
if (state[0]!=rlxd_size())
error(5);
for (k=0;k<96;k++)
{
if ((state[k+1]<0)||(state[k+1]>=167777216))
error(5);
rlxd_x.num[k]=(float)(ldexp((double)(state[k+1]),-24));
}
if (((state[97]!=0)&&(state[97]!=1))||
((state[98]!=0)&&(state[98]!=1))||
((state[99]!=0)&&(state[99]!=1))||
((state[100]!=0)&&(state[100]!=1)))
error(5);
carry.c1=(float)(ldexp((double)(state[97]),-24));
carry.c2=(float)(ldexp((double)(state[98]),-24));
carry.c3=(float)(ldexp((double)(state[99]),-24));
carry.c4=(float)(ldexp((double)(state[100]),-24));
rlxd_pr=state[101];
ir=state[102];
jr=state[103];
is=state[104];
is_old=8*ir;
prm=rlxd_pr%12;
init=1;
if (((rlxd_pr!=202)&&(rlxd_pr!=397))||
(ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
(is<0)||(is>91))
error(5);
}
#else
#define BASE 0x1000000
#define MASK 0xffffff
#define STEP(pi,pj) \
d=(*pj).c1.c1-(*pi).c1.c1-carry.c1; \
(*pi).c2.c1+=(d<0); \
d+=BASE; \
(*pi).c1.c1=d&MASK; \
d=(*pj).c1.c2-(*pi).c1.c2-carry.c2; \
(*pi).c2.c2+=(d<0); \
d+=BASE; \
(*pi).c1.c2=d&MASK; \
d=(*pj).c1.c3-(*pi).c1.c3-carry.c3; \
(*pi).c2.c3+=(d<0); \
d+=BASE; \
(*pi).c1.c3=d&MASK; \
d=(*pj).c1.c4-(*pi).c1.c4-carry.c4; \
(*pi).c2.c4+=(d<0); \
d+=BASE; \
(*pi).c1.c4=d&MASK; \
d=(*pj).c2.c1-(*pi).c2.c1; \
carry.c1=(d<0); \
d+=BASE; \
(*pi).c2.c1=d&MASK; \
d=(*pj).c2.c2-(*pi).c2.c2; \
carry.c2=(d<0); \
d+=BASE; \
(*pi).c2.c2=d&MASK; \
d=(*pj).c2.c3-(*pi).c2.c3; \
carry.c3=(d<0); \
d+=BASE; \
(*pi).c2.c3=d&MASK; \
d=(*pj).c2.c4-(*pi).c2.c4; \
carry.c4=(d<0); \
d+=BASE; \
(*pi).c2.c4=d&MASK
void RandGen::RanLxd::error(int no) const
{
switch(no)
{
case 0:
LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd");
break;
case 1:
LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)");
break;
case 2:
LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)");
break;
case 3:
LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)");
case 4:
LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd");
break;
case 5:
LATAN_ERROR(Logic, "Unexpected input data");
break;
}
}
void RandGen::RanLxd::update(void)
{
int k,kmax,d;
rlxd_dble_vec_t *pmin,*pmax,*pi,*pj;
kmax=rlxd_pr;
pmin=&rlxd_x.vec[0];
pmax=pmin+12;
pi=&rlxd_x.vec[ir];
pj=&rlxd_x.vec[jr];
for (k=0;k<kmax;k++)
{
STEP(pi,pj);
pi+=1;
pj+=1;
if (pi==pmax)
pi=pmin;
if (pj==pmax)
pj=pmin;
}
ir+=prm;
jr+=prm;
if (ir>=12)
ir-=12;
if (jr>=12)
jr-=12;
is=8*ir;
is_old=is;
}
void RandGen::RanLxd::define_constants(void)
{
int k;
one_bit=ldexp(1.0,-24);
for (k=0;k<96;k++)
{
next[k]=(k+1)%96;
if ((k%4)==3)
next[k]=(k+5)%96;
}
}
void RandGen::RanLxd::rlxd_init(int level,int seed)
{
int i,k,l;
int ibit,jbit,xbit[31];
int ix,iy;
if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
(DBL_MANT_DIG<48))
error(0);
define_constants();
if (level==1)
rlxd_pr=202;
else if (level==2)
rlxd_pr=397;
else
error(1);
i=seed;
for (k=0;k<31;k++)
{
xbit[k]=i%2;
i/=2;
}
if ((seed<=0)||(i!=0))
error(2);
ibit=0;
jbit=18;
for (i=0;i<4;i++)
{
for (k=0;k<24;k++)
{
ix=0;
for (l=0;l<24;l++)
{
iy=xbit[ibit];
ix=2*ix+iy;
xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
ibit=(ibit+1)%31;
jbit=(jbit+1)%31;
}
if ((k%4)!=i)
ix=16777215-ix;
rlxd_x.num[4*k+i]=ix;
}
}
carry.c1=0;
carry.c2=0;
carry.c3=0;
carry.c4=0;
ir=0;
jr=7;
is=91;
is_old=0;
prm=rlxd_pr%12;
init=1;
}
void RandGen::RanLxd::ranlxd(double r[],int n)
{
int k;
if (init==0)
rlxd_init(1,1);
for (k=0;k<n;k++)
{
is=next[is];
if (is==is_old)
update();
r[k]=one_bit*((double)(rlxd_x.num[is+4])+one_bit*(double)(rlxd_x.num[is]));
}
}
int RandGen::RanLxd::rlxd_size(void) const
{
return(RLXG_STATE_SIZE);
}
void RandGen::RanLxd::rlxd_get(int state[]) const
{
int k;
if (init==0)
error(3);
state[0]=rlxd_size();
for (k=0;k<96;k++)
state[k+1]=rlxd_x.num[k];
state[97]=carry.c1;
state[98]=carry.c2;
state[99]=carry.c3;
state[100]=carry.c4;
state[101]=rlxd_pr;
state[102]=ir;
state[103]=jr;
state[104]=is;
}
void RandGen::RanLxd::rlxd_reset(const int state[])
{
int k;
if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
(DBL_MANT_DIG<48))
error(4);
define_constants();
if (state[0]!=rlxd_size())
error(5);
for (k=0;k<96;k++)
{
if ((state[k+1]<0)||(state[k+1]>=167777216))
error(5);
rlxd_x.num[k]=state[k+1];
}
if (((state[97]!=0)&&(state[97]!=1))||
((state[98]!=0)&&(state[98]!=1))||
((state[99]!=0)&&(state[99]!=1))||
((state[100]!=0)&&(state[100]!=1)))
error(5);
carry.c1=state[97];
carry.c2=state[98];
carry.c3=state[99];
carry.c4=state[100];
rlxd_pr=state[101];
ir=state[102];
jr=state[103];
is=state[104];
is_old=8*ir;
prm=rlxd_pr%12;
init=1;
if (((rlxd_pr!=202)&&(rlxd_pr!=397))||
(ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
(is<0)||(is>91))
error(5);
}
#endif
// constructors ////////////////////////////////////////////////////////////////
RandGen::RandGen(void)
{
generator_.rlxd_init(RLXD_LEVEL, (int)time(NULL));
}
RandGen::RandGen(const int seed)
{
generator_.rlxd_init(RLXD_LEVEL, seed);
}
RandGen::RandGen(const RandGenState &state)
{
setState(state);
}
// state management ////////////////////////////////////////////////////////////
RandGenState RandGen::getState(void) const
{
RandGenState state;
generator_.rlxd_get(state.data());
return state;
}
void RandGen::setState(const RandGenState &state)
{
generator_.rlxd_reset(state.data());
}
// generators //////////////////////////////////////////////////////////////////
double RandGen::uniform(const double a, const double b)
{
double rx;
generator_.ranlxd(&rx, 1);
return (b-a)*rx + a;
}
unsigned int RandGen::discreteUniform(const unsigned int n)
{
return ((unsigned int)(uniform()*(double)(n)));
}
double RandGen::gaussian(const double mean, const double sigma)
{
double rx, ry, sqNrm;
do
{
rx = uniform(-1.0, 1.0);
ry = uniform(-1.0, 1.0);
sqNrm = rx*rx + ry*ry;
} while ((sqNrm > 1.0)||(sqNrm == 0.0));
return sigma*rx*sqrt(-2.0*log(sqNrm)/sqNrm) + mean;
}

103
lib/RandGen.hpp
View File

@ -1,103 +0,0 @@
/*
* RandGen.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_RandGen_hpp_
#define Latan_RandGen_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/IoObject.hpp>
#define RLXG_STATE_SIZE 105u
BEGIN_LATAN_NAMESPACE
class RandGenState: public Array<int, RLXG_STATE_SIZE, 1>,
public IoObject
{
private:
typedef Array<int, RLXG_STATE_SIZE, 1> Base;
public:
// destructor
virtual ~RandGenState(void) = default;
// IO type
IoType getType(void) const;
};
/******************************************************************************
* Random generator class *
******************************************************************************/
class RandGen
{
private:
// Martin Luescher's ranlxd generator interface
class RanLxd
{
private:
typedef struct alignas(16)
{
float c1,c2,c3,c4;
} rlxd_vec_t;
typedef struct alignas(16)
{
rlxd_vec_t c1,c2;
} rlxd_dble_vec_t;
public:
RanLxd(void);
~RanLxd(void) = default;
void ranlxd(double r[],int n);
void rlxd_init(int level,int seed);
int rlxd_size(void) const;
void rlxd_get(int state[]) const;
void rlxd_reset(const int state[]);
private:
void error(int no) const;
void update(void);
void define_constants(void);
private:
int init{0}, rlxd_pr, prm, ir, jr, is, is_old, next[96];
rlxd_vec_t one_sse, one_bit_sse, carry;
double one_bit;
union alignas(16)
{
rlxd_dble_vec_t vec[12];
float num[96];
} rlxd_x;
};
public:
// constructors
RandGen(void);
explicit RandGen(const int seed);
explicit RandGen(const RandGenState &state);
// destructor
virtual ~RandGen(void) = default;
// state management
RandGenState getState(void) const;
void setState(const RandGenState &state);
// generators
double uniform(const double a = 0.0, const double b = 1.0);
unsigned int discreteUniform(const unsigned int n);
double gaussian(const double mean = 0.0, const double sigma = 1.0);
private:
RanLxd generator_;
};
END_LATAN_NAMESPACE
#endif // Latan_RandGen_hpp_

5
utils/Makefile.am
View File

@ -7,17 +7,12 @@ endif
endif
bin_PROGRAMS = \
latan_create_rg_state \
latan_make_fake_sample\
latan_sample_combine \
latan_sample_plot_corr\
latan_sample_read \
latan_resample
latan_create_rg_state_SOURCES = create_rg_state.cpp
latan_create_rg_state_CXXFLAGS = $(COM_CXXFLAGS)
latan_create_rg_state_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_make_fake_sample_SOURCES = make_fake_sample.cpp
latan_make_fake_sample_CXXFLAGS = $(COM_CXXFLAGS)
latan_make_fake_sample_LDFLAGS = -L../lib/.libs -lLatAnalyze

44
utils/create_rg_state.cpp
View File

@ -1,44 +0,0 @@
/*
* create_rg_state.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/RandGen.hpp>
using namespace Latan;
using namespace std;
int main(int argc, char *argv[])
{
string outFilename;
if (argc != 2)
{
cerr << "usage: " << argv[0] << " <output file>" << endl;
return EXIT_FAILURE;
}
outFilename = argv[1];
RandGen gen;
Io::save(gen.getState(), outFilename);
return EXIT_SUCCESS;
}

9
utils/make_fake_sample.cpp
View File

@ -19,7 +19,6 @@
#include <iostream>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/RandGen.hpp>
using namespace std;
using namespace Latan;
@ -42,8 +41,10 @@ int main(int argc, char *argv[])
nSample = strTo<Index>(argv[3]);
outFileName = argv[4];
RandGen gen;
DSample res(nSample);
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis(val, err);
DSample res(nSample);
FOR_STAT_ARRAY(res, s)
{
@ -53,7 +54,7 @@ int main(int argc, char *argv[])
}
else
{
res[s] = gen.gaussian(val, err);
res[s] = dis(gen);
}
}
Io::save<DSample>(res, outFileName);

26
utils/resample.cpp
View File

@ -34,7 +34,7 @@ using namespace Latan;
static void usage(const string &cmdName)
{
cerr << "usage: " << cmdName;
cerr << " [-n <nsample> -b <bin size> -r <state> -o <output dir> -f {h5|sample}]";
cerr << " [-n <nsample> -b <bin size> -r <seed> -o <output dir> -f {h5|sample}]";
cerr << " <data list> <name list>";
cerr << endl;
exit(EXIT_FAILURE);
@ -43,10 +43,12 @@ static void usage(const string &cmdName)
int main(int argc, char *argv[])
{
// argument parsing ////////////////////////////////////////////////////////
int c;
string manFileName, nameFileName, stateFileName, cmdName, outDirName = ".";
string ext = "h5";
Index binSize = 1, nSample = DEF_NSAMPLE;
int c;
random_device rd;
SeedType seed = rd();
string manFileName, nameFileName, cmdName, outDirName = ".";
string ext = "h5";
Index binSize = 1, nSample = DEF_NSAMPLE;
opterr = 0;
cmdName = basename(argv[0]);
@ -64,7 +66,7 @@ int main(int argc, char *argv[])
outDirName = optarg;
break;
case 'r':
stateFileName = optarg;
seed = strTo<SeedType>(optarg);
break;
case 'f':
ext = optarg;
@ -126,25 +128,15 @@ int main(int argc, char *argv[])
// data resampling /////////////////////////////////////////////////////////
DMatSample s(nSample);
RandGen g;
RandGenState state;
cout << "-- resampling data..." << endl;
if (!stateFileName.empty())
{
state = Io::load<RandGenState>(stateFileName);
}
for (unsigned int i = 0; i < name.size(); ++i)
{
const string outFileName = name[i] + "_" + manFileName + "." + ext;
cout << '\r' << ProgressBar(i + 1, name.size());
data[name[i]].bin(binSize);
if (!stateFileName.empty())
{
g.setState(state);
}
s = data[name[i]].bootstrapMean(nSample, g);
s = data[name[i]].bootstrapMean(nSample, seed);
Io::save<DMatSample>(s, outDirName + "/" + outFileName,
File::Mode::write, outFileName);
}