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 oopse;
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:	414
Committed:	Wed Mar 9 17:58:33 2005 UTC (20 years, 7 months ago) by tim
Original Path:	*trunk/test/math/testRandNumGen.cpp*
File size:	3614 byte(s)
Log Message:	fix a bug in SectionParser (lineNo is not updated)
#	User	Rev	Content
1	tim	414	#include <functional>
2	tim	391	#include <cassert>
3			#include <fstream>
4			#include <algorithm>
5	tim	414	#include <vector>
6			#include <iostream>
7
8			#include "math/SeqRandNumGen.hpp"
9			#ifdef IS_MPI
10			#include <mpi.h>
11	tim	391	#include "math/ParallelRandNumGen.hpp"
12	tim	414	#endif
13			using namespace oopse;
14			using namespace std;
15	tim	391
16	tim	414	void testUniform(){
17			SeqRandNumGen randNumGen(823645754);
18			const int N = 100;
19	tim	391	std::vector<unsigned long int> histogram(N, 0);
20	tim	414	const unsigned long int num = 10000000;
21	tim	391	for (int i = 0; i <num; ++i) {
22	tim	414	int index = randNumGen.randInt(N -1 );
23			++histogram[index]; // rantInt returns an integer in [0, N-1]
24	tim	391	}
25	tim	414	std::ofstream uniform("uniform.dat");
26			double avg = num / N;
27	tim	391	double tolerance = 0.01*avg;
28	tim	414	for (int i = 0; i < N; ++i) {
29			//assert((histogram[i] - avg) /avg <= tolerance);
30	tim	391	uniform << i << "\t" << histogram[i] << std::endl;
31			}
32			}
33
34	tim	414	void testGaussian(){
35			SeqRandNumGen randNumGen(823645754);
36	tim	391	double mean = 100.0;
37			double variance = 1.0;
38	tim	414	const unsigned long int num = 1000000;
39	tim	391	double interval = 0.1;
40			const int size = 2000;
41			vector<unsigned long int> histogram(size , 0);
42			vector<double> normalizedHistogram(size);
43	tim	414	for (unsigned long int i = 0; i < num; ++i) {
44	tim	391	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	tim	414	for (int i = 0; i < size; ++i) {
51			gaussian << i*interval << "\t" << normalizedHistogram[i] << std::endl;
52	tim	391	}
53			}
54	tim	414	#ifdef IS_MPI
55			void testParallelRandNumGen(){
56			const unsigned long int seed = 324271632;
57			const unsigned long int nloops = 1000000;
58	tim	391	MPI_Status istatus;
59			ParallelRandNumGen mpiRandNumGen(seed);
60			const int masterNode = 0;
61	tim	414	int myRank;
62			MPI_Comm_rank( MPI_COMM_WORLD, &myRank );
63			if (myRank == masterNode) {
64	tim	391
65	tim	414	SeqRandNumGen singleRandNumGen(seed);
66			std::ofstream singleOs("single.dat");
67			std::ofstream parallelOs("parallel.dat");
68	tim	391	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	tim	414	for (unsigned long int i = 0; i < nloops; ++i) {
74	tim	391	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	tim	414	for (int j = 0; j < nProcessors; ++j) {
85			singleOs << singleRandNums[j] << "\n";
86			parallelOs << singleRandNums[j] << "\n";
87			}
88	tim	391	}
89
90
91
92			} else {
93
94			unsigned long int randNum;
95	tim	414	for (unsigned long int i = 0; i < nloops; ++i) {
96	tim	391	randNum = mpiRandNumGen.randInt();
97			MPI_Send(&randNum, 1, MPI_INT, masterNode, i, MPI_COMM_WORLD);
98			}
99
100			}
101
102			}
103	tim	414	#endif
104	tim	391	int main(int argc, char* argv[]) {
105	tim	414	#ifdef IS_MPI
106			MPI_Init(&argc, &argv);
107			std::cout << "begin test" << std::endl;
108	tim	391	if (worldRank == 0 ) {
109			testUniform();
110			testGaussian();
111			}
112
113			testParallelRandNumGen();
114
115			MPI_Finalize();
116	tim	414	#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	tim	391	}