test/math/testRandNumGen.cpp

#include <functional>
#include <cassert>
#include <fstream>
#include <algorithm>
#include <vector>
#include <iostream>

#include "math/SeqRandNumGen.hpp"
#ifdef IS_MPI
#include <mpi.h>
#include "math/ParallelRandNumGen.hpp"
#endif
using namespace OpenMD;
using namespace std;

void testUniform(){
    SeqRandNumGen randNumGen(823645754);
    const int N = 100;
    std::vector<unsigned long int> histogram(N, 0);
    const unsigned long int num = 10000000;
    for (int i = 0; i <num; ++i) {
        int index = randNumGen.randInt(N -1 );   
        ++histogram[index]; // rantInt returns an integer in [0, N-1]        
    }
    std::ofstream uniform("uniform.dat");
    double avg = num / N;
    double tolerance = 0.01*avg;
    for (int i = 0; i < N; ++i) {
        //assert((histogram[i] - avg) /avg <= tolerance);
        uniform << i << "\t" << histogram[i] << std::endl;
    }
}

void testGaussian(){      
    SeqRandNumGen randNumGen(823645754);
    double mean = 100.0;
    double variance = 1.0;
    const unsigned long int num = 1000000;
    double interval = 0.1;
    const int size = 2000;
    vector<unsigned long int> histogram(size , 0);
    vector<double> normalizedHistogram(size);
    for (unsigned long int i = 0; i < num; ++i) {
        int index = static_cast<int>(randNumGen.randNorm(mean, variance) / interval);
        ++histogram[index];        
    }

    std::transform(histogram.begin(), histogram.end(), normalizedHistogram.begin(), std::bind2nd(std::divides<double>(), num));
    std::ofstream gaussian("gaussian.dat");
    for (int i = 0; i < size; ++i) {
        gaussian << i*interval << "\t" << normalizedHistogram[i] << std::endl;
    }    
}
#ifdef IS_MPI
void testParallelRandNumGen(){      
    const unsigned long int seed = 324271632;
    const unsigned long int nloops = 1000000;
    MPI_Status istatus;
    ParallelRandNumGen mpiRandNumGen(seed);
    const int masterNode = 0;
    int myRank;
    MPI_Comm_rank( MPI_COMM_WORLD, &myRank );
    if (myRank == masterNode) {

        SeqRandNumGen singleRandNumGen(seed);
        std::ofstream singleOs("single.dat");
        std::ofstream parallelOs("parallel.dat");
        int nProcessors;
        MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
        std::vector<unsigned long int> mpiRandNums(nProcessors);
        std::vector<unsigned long int> singleRandNums(nProcessors);

        for (unsigned long int i = 0; i < nloops; ++i) {
            mpiRandNums[masterNode] = mpiRandNumGen.randInt();
        
            for (int j = 0; j < nProcessors; ++j) {
                if (j != masterNode) {
                    MPI_Recv(&mpiRandNums[j], 1, MPI_UNSIGNED_LONG, j, i, MPI_COMM_WORLD, &istatus);
                }

                singleRandNums[j] = mpiRandNumGen.randInt();
            }

            for (int j = 0; j < nProcessors; ++j) {
                singleOs << singleRandNums[j] << "\n";
                parallelOs << singleRandNums[j] << "\n";
            }
        }


    } else {

        unsigned long int randNum;
        for (unsigned long int i = 0; i < nloops; ++i) {
            randNum = mpiRandNumGen.randInt();
            MPI_Send(&randNum, 1, MPI_INT, masterNode, i, MPI_COMM_WORLD);
        }

    }

}
#endif
int main(int argc, char* argv[]) {
#ifdef IS_MPI
    MPI_Init(&argc, &argv);
    std::cout << "begin test" << std::endl;
    if (worldRank == 0 ) {
        testUniform();
        testGaussian();
    }

    testParallelRandNumGen();
    
    MPI_Finalize();
#else
   std::cout << "begin test" <<std::endl;       
   testUniform();
   testGaussian();
#endif  
   std::cout << "test done" << std::endl;  
   return 0;
}
Revision:	1465
Committed:	Fri Jul 9 23:08:25 2010 UTC (14 years, 9 months ago) by chuckv
File size:	3615 byte(s)
Log Message:	Creating busticated version of OpenMD
#	Content
1	#include <functional>
2	#include <cassert>
3	#include <fstream>
4	#include <algorithm>
5	#include <vector>
6	#include <iostream>
7
8	#include "math/SeqRandNumGen.hpp"
9	#ifdef IS_MPI
10	#include <mpi.h>
11	#include "math/ParallelRandNumGen.hpp"
12	#endif
13	using namespace OpenMD;
14	using namespace std;
15
16	void testUniform(){
17	SeqRandNumGen randNumGen(823645754);
18	const int N = 100;
19	std::vector<unsigned long int> histogram(N, 0);
20	const unsigned long int num = 10000000;
21	for (int i = 0; i <num; ++i) {
22	int index = randNumGen.randInt(N -1 );
23	++histogram[index]; // rantInt returns an integer in [0, N-1]
24	}
25	std::ofstream uniform("uniform.dat");
26	double avg = num / N;
27	double tolerance = 0.01*avg;
28	for (int i = 0; i < N; ++i) {
29	//assert((histogram[i] - avg) /avg <= tolerance);
30	uniform << i << "\t" << histogram[i] << std::endl;
31	}
32	}
33
34	void testGaussian(){
35	SeqRandNumGen randNumGen(823645754);
36	double mean = 100.0;
37	double variance = 1.0;
38	const unsigned long int num = 1000000;
39	double interval = 0.1;
40	const int size = 2000;
41	vector<unsigned long int> histogram(size , 0);
42	vector<double> normalizedHistogram(size);
43	for (unsigned long int i = 0; i < num; ++i) {
44	int index = static_cast<int>(randNumGen.randNorm(mean, variance) / interval);
45	++histogram[index];
46	}
47
48	std::transform(histogram.begin(), histogram.end(), normalizedHistogram.begin(), std::bind2nd(std::divides<double>(), num));
49	std::ofstream gaussian("gaussian.dat");
50	for (int i = 0; i < size; ++i) {
51	gaussian << i*interval << "\t" << normalizedHistogram[i] << std::endl;
52	}
53	}
54	#ifdef IS_MPI
55	void testParallelRandNumGen(){
56	const unsigned long int seed = 324271632;
57	const unsigned long int nloops = 1000000;
58	MPI_Status istatus;
59	ParallelRandNumGen mpiRandNumGen(seed);
60	const int masterNode = 0;
61	int myRank;
62	MPI_Comm_rank( MPI_COMM_WORLD, &myRank );
63	if (myRank == masterNode) {
64
65	SeqRandNumGen singleRandNumGen(seed);
66	std::ofstream singleOs("single.dat");
67	std::ofstream parallelOs("parallel.dat");
68	int nProcessors;
69	MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
70	std::vector<unsigned long int> mpiRandNums(nProcessors);
71	std::vector<unsigned long int> singleRandNums(nProcessors);
72
73	for (unsigned long int i = 0; i < nloops; ++i) {
74	mpiRandNums[masterNode] = mpiRandNumGen.randInt();
75
76	for (int j = 0; j < nProcessors; ++j) {
77	if (j != masterNode) {
78	MPI_Recv(&mpiRandNums[j], 1, MPI_UNSIGNED_LONG, j, i, MPI_COMM_WORLD, &istatus);
79	}
80
81	singleRandNums[j] = mpiRandNumGen.randInt();
82	}
83
84	for (int j = 0; j < nProcessors; ++j) {
85	singleOs << singleRandNums[j] << "\n";
86	parallelOs << singleRandNums[j] << "\n";
87	}
88	}
89
90
91
92	} else {
93
94	unsigned long int randNum;
95	for (unsigned long int i = 0; i < nloops; ++i) {
96	randNum = mpiRandNumGen.randInt();
97	MPI_Send(&randNum, 1, MPI_INT, masterNode, i, MPI_COMM_WORLD);
98	}
99
100	}
101
102	}
103	#endif
104	int main(int argc, char* argv[]) {
105	#ifdef IS_MPI
106	MPI_Init(&argc, &argv);
107	std::cout << "begin test" << std::endl;
108	if (worldRank == 0 ) {
109	testUniform();
110	testGaussian();
111	}
112
113	testParallelRandNumGen();
114
115	MPI_Finalize();
116	#else
117	std::cout << "begin test" <<std::endl;
118	testUniform();
119	testGaussian();
120	#endif
121	std::cout << "test done" << std::endl;
122	return 0;
123	}