sysrap/tests/S4MTRandGaussQTest.cc

0001 #include "S4MTRandGaussQ.h"
0002 #include "NP.hh"
0003 #include "njuffa_erfcinvf.h"
0004
0005 void test_transformQuick()
0006 {
0007     int N = 1000 ;
0008     NP* a = NP::Make<double>( N, 2 );
0009     double* aa = a->values<double>() ;
0010
0011     for(int i=0 ; i < N ; i++)
0012     {
0013         double r = double(i)/double(N-1)  ;
0014         double v = S4MTRandGaussQ::transformQuick(r) ;
0015         aa[i*2+0] = r ;
0016         aa[i*2+1] = v ;
0017     }
0018     a->save("$FOLD/test_transformQuick.npy");
0019 }
0020
0021 void test_transformQuick_vs_njuffa_erfcinvf()
0022 {
0023     int N = 1000 ;
0024
0025     int ni = N ;
0026     int nj = 3 ;
0027
0028     NP* a = NP::Make<double>( ni, nj );
0029     double* aa = a->values<double>() ;
0030
0031     for(int i=0 ; i < ni ; i++)
0032     {
0033         double r = double(i)/double(N-1)  ;
0034         aa[i*nj+0] = r ;
0035         aa[i*nj+1] = S4MTRandGaussQ::transformQuick(r)  ;
0036         aa[i*nj+2] = -sqrtf(2.f)*njuffa_erfcinvf(r*2.f) ;
0037         // argument to erfcinvf needs to be in range 0->2
0038         // then have to scale to match S4MTRandGaussQ::transformQuick
0039     }
0040     a->save("$FOLD/test_transformQuick_vs_njuffa_erfcinvf.npy");
0041 }
0042
0043
0044 void test_shoot()
0045 {
0046     int N = 1000000 ;
0047     NP* a = NP::Make<double>( N );
0048     double* aa = a->values<double>() ;
0049
0050     double mean = 0.0 ;
0051     double stdDev = 0.1 ;
0052
0053     for(int i=0 ; i < N ; i++) aa[i] = S4MTRandGaussQ::shoot(mean, stdDev) ;
0054     a->save("$FOLD/test_shoot.npy");
0055 }
0056
0057
0058 int main()
0059 {
0060     //test_transformQuick();
0061     //test_shoot();
0062
0063     test_transformQuick_vs_njuffa_erfcinvf();
0064
0065     return 0 ;
0066 }