// program description
const char programName[21] = "genePed";
    // program name
const char versionNumber[15] = "2.0";
    // version number
const char versionDate[15] = "11/30/2006";
    // version date
const char author[61] = "Chih-Horng Kuo and John Avise";
    // author
const char email[21] = "chkuo@lifedev.org";
    // author email
const char web[61] = "http://chkuo.name/";
    // author web

// This program is developed to simulate 
// gene trees from organismal pedigree

// preprocessor directive
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <string>
#include <ctype.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

// For standard Template Library (STL)
using namespace std;

// constants 
    // filename length
	const int fileNameLength = 81;
    // id string length
	const int idStringLength = 21;
    // comment line length
	const int commentLength = 121;
    // maxRandomInt
    // controls precision
	const int maxRandomInt = 10000;
    // number of digits after the decimal point
    // used to control display/output file    
	const int nDigit = 4;
    // report interval    
	const int repInterval = 50;
    // symbol in input/output file for comment    
	const char symbolComment = '/';
    // symbol in input/output file for nPopSize    
	const char symbolNPopSize = '#';
    // symbol in input/output file for nDispStep    
	const char symbolNDispStep = '^';
    // symbol in input/output file for the start of dispersal prob table
	const char symbolDispProb = '@';
    
    
// typedef
    // new data type "fileName" to hold file name
	typedef char fileName[fileNameLength];
    // new data type "idString"  to hold id strings
	typedef char idString[idStringLength];
    
// start of program functions

// welcome function
// display welcome info
void welcome(void) {
    cout << "Program name:\t" << programName << endl;
    cout << "\tVersion number:\t" << versionNumber << endl;
    cout << "\tVersion date:\t" << versionDate << endl;
    cout << endl;

    cout << "Contact information" << endl; 
    cout << "\tAuthor:\t" << author << endl; 
    cout << "\temail:\t" << email << endl;
    cout << "\tWeb:\t" << web << endl;
    cout << endl;
} // end of welcome function

// memoryError function
void memoryError(void) {
    cout << endl;
    cout << "*** Error allocating memory!" << endl;
    exit(0);                     
} // end of memoryError function

// openInputFileDsp function
// open input file and get input file name
void openInputFileDsp(ifstream &inputFileDsp, fileName &inputNameDsp) {
 
    // prompt for input filename
    cout << endl;
    cout << "Input file must be in the same directory as this program," << endl;
    cout << "The limitation of file name length is " << (fileNameLength - 1) << " characters." << endl;
    cout << "Please enter file name of the input file (Dispersal prob <.dsp>): ";
    cin >> inputNameDsp;

	// open input file
	inputFileDsp.open(inputNameDsp);
	// check input file
	if ( inputFileDsp == 0 ) {
	    cout << endl;
		cout << "*** Error opening input file." << endl;
		exit(0);
	}
	else
		cout << endl;
		cout << "Input file (" << inputNameDsp << ") opened successfully." << endl;

} // end of openInputFileDsp


// openOutputRep function
// open output file and get output file name
void openOutputRep(ofstream &outputFileRep, fileName &outputNameRep) {
    // prompt for output filename
    cout << endl;
    cout << "The limitation of file name length is " 
         << (fileNameLength - 1) 
         << " characters."
         << endl; 
    cout << "Please enter file name of the output file (Report <.rep>): ";
    cin >> outputNameRep;
    
    // open outputFile
        outputFileRep.open(outputNameRep);
        if ( outputFileRep == 0 ) {
    	    cout << endl;
    		cout << "*** Error opening output file." << endl;
            exit(0);
        }
        else
		cout << endl;
		cout << "Output file (" << outputNameRep << ") opened successfully." << endl;

} // end of openOutputRep

// getRandomInt function
// receive minimum and maximum of the random int desired
// return random number within the min/max range
int getRandomInt(int min, int max) {
    return ((rand() % (max + 1 - min)) + min);
} // end of getRandomInt function

// getModule function
void getModule(int &module) {
    // display available modules 
    cout << endl;
    cout << "Modules available:" << endl;
    cout << "  (0) End program" << endl;
    cout << "  (1) Module 1 (4 alleles, equal spacing at [0, 0.25, 0.5, 0.75])" << endl;
    cout << "  (2) Module 2 (4 alleles, unequal spacing at [0, 0.125, 0.25, 0.5])" << endl;
    cout << "  (3) Module 3 (5 alleles, equal spacing at [0, 0.2, 0.4, 0.6, 0.8])" << endl;
    
    do {
        cout << endl
             << "Please select the module you want to use: ";
        cin >> module;
    } while (module < 0 || module > 3);
    
    cout << endl; 
    cout << "Module selected = " << module << endl;

} // getModule function

// getNIter
void getNIter(int &nIter) {
    do {
        cout << endl
             << "Please enter the number of iterations: ";
        cin >> nIter;
    } while (nIter < 1);

} // getNIter

// getNPedigree
void getNPedigree(int &nPedigree) {
    do {
        cout << endl
             << "Please enter the number of independent pedigrees: ";
        cin >> nPedigree;
    } while (nPedigree < 1);
} // getNPedigree

// getNTree
void getNTree(int &nTree) {
    do {
        cout << endl
             << "Please enter the number of gene trees per pedigree: ";
        cin >> nTree;
    } while (nTree < 1);
} // getNTree

// getDispTable
void getDispTable(int nDispStep, int *dispStep, float *dispProb, int *dispTable) {
    //
    int currentNStep;
    int flagFoundD;
    
    //
    int countD;
    int countR;
    
    // 
    for (countD = 0; countD < nDispStep; countD++) {
        dispProb[countD] = dispProb[countD] * maxRandomInt;
    }
    
    // reset
    countD = 0;
    currentNStep = dispStep[0];

    //
    for (countR = 0; countR < (maxRandomInt + 1); countR++) {
        if (countR <= dispProb[countD]) {
            dispTable[countR] = dispStep[countD];
        }
        else {
            do {
                // reset flag
                flagFoundD = 0;
                countD++;
               
                // check if index within range
                if (countD >= nDispStep) {
                    // if outside the range
                    cout << endl;
                    cout << "*** Dsp setting error ***" << endl;
                    cout << "*** index out of range ***" << endl;
                    exit(0);                     
                }
                else {
                    if (countR <= dispProb[countD]) {
                        dispTable[countR] = dispStep[countD];
                        flagFoundD = 1;
                    }
                }
            } while(flagFoundD != 1);
        }
    }
    
    
} // getDispTable

// getParLoc
int getParLoc(int selfLoc, int nPopSize, int *dispTable) {
    //
    int randInt;    // random int between 0  and maxRandomInt
    int nStepMoved; // number of steps moved
    int parLoc;     // location of the parent
    
    randInt = getRandomInt(0, maxRandomInt);
    nStepMoved = dispTable[randInt];
    
    parLoc = selfLoc + nStepMoved;
    
    // check if outside of the boundary (in a ring shaped habitat)
    // left boundary = 0, right boundary = (nPopSize - 1)
    while (parLoc < 0 || parLoc > (nPopSize - 1)) {
        if (parLoc < 0) {
            // if out of left boundary
            parLoc = nPopSize + parLoc;
        }
        else {
            // if out of right boundary
            parLoc = parLoc - nPopSize;
        }
    } // while outside of the boundary (in a ring shaped habitat)
    
    return parLoc;

} // getParLoc

// getTreeTopology4
// for trees of 4 alleles
void getTreeTopology4(int nPedigree, int nTree, int **event1, int **event2, int **treeTopology) {
    // the tree topology was determined based on 2 coal events
    // trees of 4 alleles (alleleID = 0, 1, 2, 3)
    // event and treeTopology ID table
    // event1   event2  treeTopology    Tree
    // 01       02      0               (((0,1),2),3)
    // 01       03      1               (((0,1),3),2)
    // 01       23      2               ((0,1),(2,3))
    // 02       01      3               (((0,2),1),3)
    // 02       03      4               (((0,2),3),1)
    // 02       13      5               ((0,2),(1,3))
    // 03       01      6               (((0,3),1),2)
    // 03       02      7               (((0,3),2),1)
    // 03       12      8               ((0,3),(1,2))
    // 12       01      9               (((1,2),0),3)
    // 12       13      10              (((1,2),3),0)
    // 12       03      11              ((1,2),(0,3))
    // 13       01      12              (((1,3),0),2)
    // 13       12      13              (((1,3),2),0)
    // 13       02      14              ((1,3),(0,2))
    // 23       02      15              (((2,3),0),1)
    // 23       12      16              (((2,3),1),0)
    // 23       01      17              ((2,3),(0,1))

    int countP;
    int countT;

    // get treeTopology (based on event1, event2)
    for (countP = 0; countP < nPedigree; countP++) {
        for (countT = 0; countT < nTree; countT++) {
            if (event1[countP][countT] == 1) {
                if (event2[countP][countT] == 2) {
                    treeTopology[countP][countT] = 0;
                }
                else if (event2[countP][countT] == 3) {
                    treeTopology[countP][countT] = 1;
                }
                else if (event2[countP][countT] == 23) {
                    treeTopology[countP][countT] = 2;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else if (event1[countP][countT] == 2) {
                if (event2[countP][countT] == 1) {
                    treeTopology[countP][countT] = 3;
                }
                else if (event2[countP][countT] == 3) {
                    treeTopology[countP][countT] = 4;
                }
                else if (event2[countP][countT] == 13) {
                    treeTopology[countP][countT] = 5;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else if (event1[countP][countT] == 3) {
                if (event2[countP][countT] == 1) {
                    treeTopology[countP][countT] = 6;
                }
                else if (event2[countP][countT] == 2) {
                    treeTopology[countP][countT] = 7;
                }
                else if (event2[countP][countT] == 12) {
                    treeTopology[countP][countT] = 8;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else if (event1[countP][countT] == 12) {
                if (event2[countP][countT] == 1) {
                    treeTopology[countP][countT] = 9;
                }
                else if (event2[countP][countT] == 13) {
                    treeTopology[countP][countT] = 10;
                }
                else if (event2[countP][countT] == 3) {
                    treeTopology[countP][countT] = 11;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else if (event1[countP][countT] == 13) {
                if (event2[countP][countT] == 1) {
                    treeTopology[countP][countT] = 12;
                }
                else if (event2[countP][countT] == 12) {
                    treeTopology[countP][countT] = 13;
                }
                else if (event2[countP][countT] == 2) {
                    treeTopology[countP][countT] = 14;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else if (event1[countP][countT] == 23) {
                if (event2[countP][countT] == 2) {
                    treeTopology[countP][countT] = 15;
                }
                else if (event2[countP][countT] == 12) {
                    treeTopology[countP][countT] = 16;
                }
                else if (event2[countP][countT] == 1) {
                    treeTopology[countP][countT] = 17;
                }
                else {
                    // error occured in event2!
                    cout << "event2 error!" 
                         << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                         << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
                         << endl;
                    exit(0);
                }
            }
            else {
                // error occured in event1!
            	cout << "event1 error!" 
                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                     << endl;
                exit(0);
            }
        }
    }// get treeTopology (based on event1, event2)

} // getTreeTopology4

// getTreeTopology5
// for trees of 5 alleles
void getTreeTopology5(int nPedigree, int nTree, int **event1, int **event2, int **event3, int **treeTopology) {
    // the tree topology was determined based on 3 coal events
    // trees of 5 alleles (alleleID = 0, 1, 2, 3, 4)
    // event and treeTopology ID table

/*
event1	event2	event3	treeId
01	02	03	0
01	02	04	1
01	02	34	2
01	03	02	3
01	03	04	4
01	03	24	5
01	04	02	6
01	04	03	7
01	04	34	8
01	23	02	9
01	23	04	10
01	23	24	11
01	24	02	12
01	24	03	13
01	24	23	14
01	34	02	15
01	34	03	16
01	34	23	17
02	01	03	18
02	01	04	19
02	01	34	20
02	03	01	21
02	03	04	22
02	03	14	23
02	04	01	24
02	04	03	25
02	04	13	26
02	13	01	27
02	13	04	28
02	13	14	29
02	14	01	30
02	14	03	31
02	14	13	32
02	34	01	33
02	34	03	34
02	34	13	35
03	01	02	36
03	01	04	37
03	01	24	38
03	02	01	39
03	02	04	40
03	02	14	41
03	04	01	42
03	04	02	43
03	04	12	44
03	12	01	45
03	12	04	46
03	12	14	47
03	14	01	48
03	14	02	49
03	14	12	50
03	24	01	51
03	24	02	52
03	24	12	53
04	01	02	54
04	01	03	55
04	01	23	56
04	02	01	57
04	02	03	58
04	02	13	59
04	03	01	60
04	03	02	61
04	03	12	62
04	12	01	63
04	12	03	64
04	12	13	65
04	13	01	66
04	13	02	67
04	13	12	68
04	23	01	69
04	23	02	70
04	23	12	71
12	01	03	72
12	01	04	73
12	01	34	74
12	03	01	75
12	03	04	76
12	03	14	77
12	04	01	78
12	04	03	79
12	04	13	80
12	13	01	81
12	13	04	82
12	13	14	83
12	14	01	84
12	14	03	85
12	14	13	86
12	34	01	87
12	34	03	88
12	34	13	89
13	01	02	90
13	01	04	91
13	01	24	92
13	02	01	93
13	02	04	94
13	02	14	95
13	04	01	96
13	04	02	97
13	04	12	98
13	12	01	99
13	12	04	100
13	12	14	101
13	14	01	102
13	14	02	103
13	14	12	104
13	24	01	105
13	24	02	106
13	24	12	107
14	01	02	108
14	01	03	109
14	01	23	110
14	02	01	111
14	02	03	112
14	02	13	113
14	03	01	114
14	03	02	115
14	03	12	116
14	12	01	117
14	12	03	118
14	12	13	119
14	13	01	120
14	13	02	121
14	13	12	122
14	23	01	123
14	23	02	124
14	23	12	125
23	01	02	126
23	01	04	127
23	01	24	128
23	02	01	129
23	02	04	130
23	02	14	131
23	04	01	132
23	04	02	133
23	04	12	134
23	12	01	135
23	12	04	136
23	12	14	137
23	14	01	138
23	14	02	139
23	14	12	140
23	24	01	141
23	24	02	142
23	24	12	143
24	01	02	144
24	01	03	145
24	01	23	146
24	02	01	147
24	02	03	148
24	02	13	149
24	03	01	150
24	03	02	151
24	03	12	152
24	12	01	153
24	12	03	154
24	12	13	155
24	13	01	156
24	13	02	157
24	13	12	158
24	23	01	159
24	23	02	160
24	23	12	161
34	01	02	162
34	01	03	163
34	01	23	164
34	02	01	165
34	02	03	166
34	02	13	167
34	03	01	168
34	03	02	169
34	03	12	170
34	12	01	171
34	12	03	172
34	12	13	173
34	13	01	174
34	13	02	175
34	13	12	176
34	23	01	177
34	23	02	178
34	23	12	179
*/

    int countP;
    int countT;

    // get treeTopology (based on event1, event2, event3)
    for (countP = 0; countP < nPedigree; countP++) {
        for (countT = 0; countT < nTree; countT++) {
            if (event1[countP][countT] == 1) {
	            if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 0;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 1;
		            }
		            else if (event3[countP][countT] == 34) {
	                    treeTopology[countP][countT] = 2;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 3;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 4;
		            }
		            else if (event3[countP][countT] == 24) {
	                    treeTopology[countP][countT] = 5;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 6;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 7;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 8;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 23) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 9;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 10;
		            }
		            else if (event3[countP][countT] == 24) {
	                    treeTopology[countP][countT] = 11;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 24) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 12;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 13;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 14;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 34) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 15;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 16;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 17;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 2) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 18;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 19;
		            }
		            else if (event3[countP][countT] == 34) {
	                    treeTopology[countP][countT] = 20;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 21;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 22;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 23;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 24;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 25;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 26;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 27;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 28;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 29;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 14) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 30;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 31;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 32;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 34) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 33;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 34;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 35;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 3) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 36;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 37;
		            }
		            else if (event3[countP][countT] == 24) {
	                    treeTopology[countP][countT] = 38;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 39;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 40;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 41;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 42;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 43;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 44;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 45;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 46;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 47;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 14) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 48;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 49;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 50;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 24) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 51;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 52;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 53;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 4) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 54;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 55;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 56;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 57;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 58;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 59;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 60;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 61;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 62;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 63;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 64;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 65;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 66;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 67;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 68;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 23) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 69;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 70;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 71;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 12) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 72;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 73;
		            }
		            else if (event3[countP][countT] == 34) {
	                    treeTopology[countP][countT] = 74;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 75;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 76;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 77;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 78;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 79;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 80;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 81;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 82;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 83;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 14) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 84;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 85;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 86;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 34) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 87;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 88;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 89;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 13) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 90;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 91;
		            }
		            else if (event3[countP][countT] == 24) {
	                    treeTopology[countP][countT] = 92;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 93;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 94;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 95;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 96;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 97;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 98;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 99;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 100;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 101;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 14) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 102;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 103;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 104;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 24) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 105;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 106;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 107;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 14) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 108;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 109;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 110;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 111;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 112;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 113;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 114;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 115;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 116;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 117;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 118;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 119;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 120;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 121;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 122;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 23) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 123;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 124;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 125;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 23) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 126;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 127;
		            }
		            else if (event3[countP][countT] == 24) {
	                    treeTopology[countP][countT] = 128;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 129;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 130;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 131;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 4) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 132;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 133;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 134;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 135;
		            }
		            else if (event3[countP][countT] == 4) {
	                    treeTopology[countP][countT] = 136;
		            }
		            else if (event3[countP][countT] == 14) {
	                    treeTopology[countP][countT] = 137;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 14) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 138;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 139;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 140;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 24) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 141;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 142;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 143;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 24) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 144;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 145;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 146;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 147;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 148;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 149;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 150;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 151;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 152;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 153;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 154;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 155;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 156;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 157;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 158;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 23) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 159;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 160;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 161;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else if (event1[countP][countT] == 34) {
	            if (event2[countP][countT] == 1) {
		            if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 162;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 163;
		            }
		            else if (event3[countP][countT] == 23) {
	                    treeTopology[countP][countT] = 164;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 2) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 165;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 166;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 167;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 3) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 168;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 169;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 170;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 12) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 171;
		            }
		            else if (event3[countP][countT] == 3) {
	                    treeTopology[countP][countT] = 172;
		            }
		            else if (event3[countP][countT] == 13) {
	                    treeTopology[countP][countT] = 173;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 13) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 174;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 175;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 176;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else if (event2[countP][countT] == 23) {
		            if (event3[countP][countT] == 1) {
	                    treeTopology[countP][countT] = 177;
		            }
		            else if (event3[countP][countT] == 2) {
	                    treeTopology[countP][countT] = 178;
		            }
		            else if (event3[countP][countT] == 12) {
	                    treeTopology[countP][countT] = 179;
		            }
		            else {
		                // error occured in event3!
		            	cout << "event3 error!" 
		                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
		                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
		                     << "\tevent3[" << countP << "][" << countT << "] = "<< event3[countP][countT] 
		                     << endl;
		                exit(0);
		            }
	            }
	            else {
	                // error occured in event2!
	            	cout << "event2 error!" 
	                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
	                     << "\tevent2[" << countP << "][" << countT << "] = "<< event2[countP][countT] 
	                     << endl;
	                exit(0);
	            }
            }
            else {
                // error occured in event1!
            	cout << "event1 error!" 
                     << "\tevent1[" << countP << "][" << countT << "] = "<< event1[countP][countT] 
                     << endl;
                exit(0);
            }
        }
    }// get treeTopology (based on event1, event2)

} // getTreeTopology5


// statTopology
void statTopology(int countI, int nTopology, int nPedigree, int nTree, int **treeTopology, int **topologyCount, float **topologyFreq) {
    int countP;
    int countT;
    
    // reset
    for (countT = 0; countT < nTopology; countT++) {
        topologyCount[countI][countT] = 0;
    }
    
    // get count for each tree
    for (countP = 0; countP < nPedigree; countP++) {
        for (countT = 0; countT < nTree; countT++) {
            if (treeTopology[countP][countT] < 0 || treeTopology[countP][countT] >= nTopology) {
                cout << "treeTopology id error!\t";
                cout << "treeTopology[" << countP << "][" << countT << "] = " << treeTopology[countP][countT] << endl;
                exit(0);
            }
            else {
                topologyCount[countI][(treeTopology[countP][countT])]++;
            }
        }
    }    
    
    // get freq
    for (countT = 0; countT < nTopology; countT++) {
        topologyFreq[countI][countT] = ( float (topologyCount[countI][countT]) / (nPedigree * nTree));
        //cout << "topologyCount[" << countI << "][" << countT << "] = " << topologyCount[countI][countT] << endl;
        //cout << "topologyFreq[" << countI << "][" << countT << "] = " << topologyFreq[countI][countT] << endl;
    }

} // statTopology

// statTopologyFreq
void statTopologyFreq(int countI, int nTopology, float **topologyFreq, float *topologyFreqMax, float *topologyFreqMin, float *topologyFreqVar)
{
    int countT;
    float sum = 0;

    // find max/min
    // reset 
    topologyFreqMax[countI] = topologyFreq[countI][0];
    topologyFreqMin[countI] = topologyFreq[countI][0];
    for (countT = 0; countT < nTopology; countT++) {
		sum += topologyFreq[countI][countT];
        // update max
        if (topologyFreqMax[countI] < topologyFreq[countI][countT])
            topologyFreqMax[countI] = topologyFreq[countI][countT];
        // update min
        if (topologyFreqMin[countI] > topologyFreq[countI][countT])
            topologyFreqMin[countI] = topologyFreq[countI][countT];
    }
    
    // get variance
    // use population variance (divide by N instead of (N-1))
    // because all possible topologies are considered here
    float mean = (sum / nTopology);
    float sumSquare = 0;
    
    for (countT = 0; countT < nTopology; countT++) {
    	sumSquare += (topologyFreq[countI][countT] - mean) * (topologyFreq[countI][countT] - mean);
    }
    topologyFreqVar[countI] = (sumSquare / nTopology);
    
} // statTopologyFreq

// statFloat
void statFloat(int nSize, float *data, float &mean, float &variance, float &stdDev, float &max, float &min) {
    // 
    float sum = 0;
    float sumSquare = 0;

    int countN;
    
    max = data[0];
    min = data[0];
    for (countN = 0; countN < nSize; countN++)
    {
        sum += data[countN];
        if (max < data[countN])
            max = data[countN];
        if (min > data[countN])
            min = data[countN];
    }
    
    mean = sum / nSize;
    
    for (countN = 0; countN < nSize; countN++) {
        sumSquare += (data[countN] - mean) * (data[countN] - mean);
    }
    variance = (sumSquare / (nSize - 1));
    stdDev = (sqrt (sumSquare / (nSize - 1)));
} // statFloat

// readInputSettingDsp function
void readInputSettingDsp(ifstream &inputFileDsp, int &nPopSize, int &nDispStep) {
    // declare variables
    char symbolInput;
    char inputComment[81];
    int flagNPopSize = 0;
    int flagNDispStep = 0;
    
    // read input setting
    do {
        inputFileDsp.get(symbolInput);
        switch (symbolInput) {
            case symbolComment:     inputFileDsp.getline(inputComment, 81, '\n');
                                    break;
            case symbolNPopSize:    inputFileDsp >> nPopSize;
                                    flagNPopSize = 1;
                                    break;
            case symbolNDispStep:   inputFileDsp >> nDispStep;
                                    flagNDispStep = 1;
                                    break;
        }
    } while (flagNPopSize !=1 || flagNDispStep !=1);

    // check input setting

    // nPopSize
    if (nPopSize < 1) {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** nPopSize < 1" << endl;
		exit(0);
    }
    if (nPopSize % 4 != 0) {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** nPopSize % 4 != 0" << endl;
		exit(0);
    }

    // nDispStep
    if (nDispStep < 1) {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** nDispStep < 1" << endl;
		exit(0);
    }

} // readInputSettingDsp function

// readInputDsp
void readInputDsp(ifstream &inputFileDsp, int nDispStep, int *dispStep, float *dispProb) {
    //
    int countD;
    // declare variables
    char symbolInput;
    char inputComment[81];
    int flagDispProb = 0;
    
    // read input setting
    do {
        inputFileDsp.get(symbolInput);
        switch (symbolInput) {
            case symbolComment:     inputFileDsp.getline(inputComment, 81, '\n');
                                    break;
            case symbolDispProb:    
                                    for (countD = 0; countD < nDispStep; countD++) {
                                        inputFileDsp >> dispStep[countD];
                                        inputFileDsp >> dispProb[countD];        
                                    }
                                    flagDispProb = 1;
                                    break;
        }
    } while (flagDispProb !=1);

    // check input
    // normalize cum prob
	for (countD = 0; countD < nDispStep; countD++) {
        dispProb[countD] = dispProb[countD] / dispProb[(nDispStep - 1)];
    }
} // readInputDsp

// dispInputSettingDsp
void dispInputSettingDsp(fileName &inputNameDsp, int nPopSize, int nDispStep, int *dispStep, float *dispProb) {
	//
	int countD;
	
	//
	cout << endl;
	cout << "Input setting from " << inputNameDsp << endl;
	cout << endl;
	
	cout << "  nPopSize = " << nPopSize << endl;
	cout << "  nDispStep = " << nDispStep << endl;
	cout << endl;
	
	cout << "  dispStep" << "\t" << "dispProb" << endl;
	for (countD = 0; countD < nDispStep; countD++) {
		cout << "  " << dispStep[countD];
		cout << "\t" << setprecision(nDigit) 
             << setiosflags(ios::fixed | ios::showpoint | ios::right)
             << dispProb[countD] << endl;
	}
	cout << endl;	

} // dispInputSettingDsp

// writeHeader function
void writeHeader() {   
    // get time
    time_t rawtime;
    struct tm * timeinfo;

    time ( &rawtime );
    timeinfo = localtime ( &rawtime );

    // write header 
    cout << "***\tStart of header\t***" << endl;
    cout << "Program name:\t" << programName << endl;
    cout << "\tVersion number:\t" << versionNumber << endl;
    cout << "\tRelease date:\t" << versionDate << endl;
    cout << endl;

    cout << "Contact information" << endl; 
    cout << "\tAuthor:\t" << author << endl; 
    cout << "\temail:\t" << email << endl;
    cout << "\tWeb:\t" << web << endl;
    cout << endl;
    
    cout << "Run date and time:\t" << asctime(timeinfo);
    cout << "***\tEnd of header\t***" << endl;

} // end of writeHeader function

// writeSetting
void writeSetting(fileName inputNameDsp, fileName outputNameRep, int nPopSize, int nDispStep, int nIter, int nPedigree, int nTree) {
    cout << endl;
    cout << "***\tStart of Settings\t***" << endl;

    cout << "Input and output files" << endl;
    cout << "\tDispersal prob file:\t" << inputNameDsp << endl;
    cout << "\tReport file:\t" << outputNameRep << endl;
    cout << endl;


    cout << "User settings" << endl; 
    cout << "\tnPopSize =\t" << nPopSize << endl;
    cout << "\tnDispStep =\t" << nDispStep << endl;
    cout << "\tnIter =\t" << nIter << endl;
    cout << "\tnPedigree =\t" << nPedigree << endl;
    cout << "\tnTree =\t" << nTree << endl;

    cout << "***\tEnd of Settings\t***" << endl;

} // writeSetting

// writeSummary
void writeSummary(int nIter, float *topologyFreqMax, float *topologyFreqMin, float *topologyFreqVar) {
    float mean;
    float variance;
    float stdDev;
    float max;
    float min;

    cout << endl;
    cout << "***\tStart of Summary\t***" << endl;
    
    // header
    cout << "Stat\tMean\tVariance\tStdDev\tMax\tMin" << endl;
    
    // topologyFreqMax
    statFloat(nIter, topologyFreqMax, mean, variance, stdDev, max, min);
    
    cout << "topologyFreqMax"
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< mean
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< variance
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< stdDev
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< max
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< min
		<< endl;
    
    // topologyFreqMin
    statFloat(nIter, topologyFreqMin, mean, variance, stdDev, max, min);
    
    cout << "topologyFreqMin"
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< mean
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< variance
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< stdDev
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< max
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< min
		<< endl;
    
    // topologyFreqMin
    statFloat(nIter, topologyFreqVar, mean, variance, stdDev, max, min);
    
    cout << "topologyFreqVar"
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< mean
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< variance
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< stdDev
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< max
		<< "\t" << setprecision(nDigit) 
		<< setiosflags(ios::fixed | ios::showpoint | ios::right)
		<< min
		<< endl;
    
    cout << "***\tEnd of Summary\t***" << endl;
} // writeSummary

// writeTopologyFreq
void writeTopologyFreq(int nIter, int nTopology, float **topologyFreq, float *topologyFreqMax, float *topologyFreqMin, float *topologyFreqVar) {
    int countI;
    int countT;
    
    cout << endl;
    cout << "***\tStart of Topology Freq\t***" << endl;

    // table header
    cout << "Iter\tMax\tMin\tVar";
    for (countT = 0; countT < nTopology; countT++) {
        cout << "\tTopo" << (countT+1);
    }
    cout << endl;
    
    // table body
    for (countI = 0; countI < nIter; countI++) {
        cout << (countI + 1);
        cout << "\t" << setprecision(nDigit) 
                   << setiosflags(ios::fixed | ios::showpoint | ios::right)
                   << topologyFreqMax[countI];
        cout << "\t" << setprecision(nDigit) 
                   << setiosflags(ios::fixed | ios::showpoint | ios::right)
                   << topologyFreqMin[countI];
        cout << "\t" << setprecision(nDigit) 
                   << setiosflags(ios::fixed | ios::showpoint | ios::right)
                   << topologyFreqVar[countI];
        for (countT = 0; countT < nTopology; countT++) {
            cout << "\t" << setprecision(nDigit) 
                       << setiosflags(ios::fixed | ios::showpoint | ios::right)
                       << topologyFreq[countI][countT];
        }
        cout << endl;
    }

    cout << "***\tEnd of Topology Freq\t***" << endl;

} // writeTopologyFreq


// writeTopologyVar
void writeTopologyVar(ofstream &outputFile, int nIter, float *topologyFreqVar) {
    int countI;
    
    for (countI = 0; countI < nIter; countI++) {
        outputFile << setprecision(nDigit)
        		   << setiosflags(ios::scientific)
                   << topologyFreqVar[countI]
                   << endl;
    }
} // writeTopologyVar


// module1
// (4 alleles, equal spacing at [0, 0.25, 0.5, 0.75]
void module1(void) {
    // nAllele
	const int nAllele = 4;
    // nTopology
	const int nTopology = 18;

    // file operation
        // input file - dispersal prob file
            // declare fileName
            fileName inputNameDsp; 
            // declare input file
            ifstream inputFileDsp;
            // call openInputFile function
            openInputFileDsp(inputFileDsp, inputNameDsp);

        // report file 
            // declare fileName
            fileName outputNameRep; 
            // declare output file
            ofstream outputFileRep;
            // call openOutput function
            openOutputRep(outputFileRep, outputNameRep);

        //
        writeHeader();

    // declare variables 
        // settings
        int nPopSize;   // pop size
        int nDispStep;  // number of possible dispersal steps 

        int nIter;      // number of iterations
        int nPedigree;  // number of pedigrees
        int nTree;      // number of gene trees per pedigree

        // variable
        int distAllele; // distance between alleles (= nPopSize / 4)
        int flagRun;    // bool flag for while loop

        // loop counting
        int countI;     // count iteration
        int countG;     // count generation
        int countN;     // count individual
        int countP;     // count pedigree
        int countT;     // count tree
        int countA;     // count allele
        int countA1;    // count allele
        int countA2;    // count allele

    // get settings from input file
    // readInputSettingDsp
    readInputSettingDsp(inputFileDsp, nPopSize, nDispStep);

    // distAllele (set to equal dist)
    distAllele = nPopSize / nAllele;
            
    // get settings from user
    getNIter(nIter);
    getNPedigree(nPedigree);
    getNTree(nTree);
    
    // write settings to stdout
    writeSetting(inputNameDsp, outputNameRep, nPopSize, nDispStep, nIter, nPedigree, nTree);
    
    // declare and initiate array   
    cout << "initiating data arrays..." << endl;     
        // array to hold dispStep
        int *dispStep;
        dispStep = new int [nDispStep];
        if (dispStep == NULL)
            memoryError();        
        
        // array to hold dispProb 
        float *dispProb;
        dispProb = new float [nDispStep];
        if (dispProb == NULL)
            memoryError();        

        // array to hold dispTable
        // hush table to find nStep
        int *dispTable;
        dispTable = new int [(maxRandomInt + 1)];
        if (dispTable == NULL)
            memoryError();
        
        // parent location
        int *p0Loc;
        p0Loc = new int [nPopSize];
        if (p0Loc == NULL)
            memoryError();
        
        // parent location
        int *p1Loc;
        p1Loc = new int [nPopSize];
        if (p1Loc == NULL)
            memoryError();
                    
        // array to hold location in current generation [nTree][nAllele]
        int **curLoc;
        curLoc = new int *[nTree];
        if (curLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			curLoc[countT] = new int [nAllele];
	        if (curLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold location in previous generation [nTree][nAllele]
        int **preLoc;
        preLoc = new int *[nTree];
        if (preLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			preLoc[countT] = new int [nAllele];
	        if (preLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold allele that has been coaled [nTree][nAllele]
        // either itself (have not coaled with any other allele)
        // or a smaller alleleId
        // when more than 2 alleles coaled, use the smallest
        //   e.g if allele 1, 2, 4 coaled
        //   allele[countT][0] = 0
        //   allele[countT][1] = 0
        //   allele[countT][3] = 0
        int **alleleCoal;
        alleleCoal = new int *[nTree];
        if (alleleCoal == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			alleleCoal[countT] = new int [nAllele];
	        if (alleleCoal == NULL)
	            memoryError();			 
		}
        
        // array to hold coal event1 id [nPedigree][nTree]
        int **event1;
        event1 = new int *[nPedigree];
        if (event1 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event1[countP] = new int [nTree];
	        if (event1 == NULL)
	            memoryError();			 
		}

        // array to hold coal event2 id [nPedigree][nTree]
        int **event2;
        event2 = new int *[nPedigree];
        if (event2 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event2[countP] = new int [nTree];
	        if (event2 == NULL)
	            memoryError();			 
		}

        // array to hold tree topology id [nPedigree][nTree]
        int **treeTopology;
        treeTopology = new int *[nPedigree];
        if (treeTopology == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			treeTopology[countP] = new int [nTree];
	        if (treeTopology == NULL)
	            memoryError();			 
		}
        
        // data summary arrays
        
        // array to hold tree topology count [nIter][nTopology]
        int **topologyCount;
        topologyCount = new int *[nIter];
        if (topologyCount == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyCount[countI] = new int [nTopology];
	        if (topologyCount == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq [nIter][nTopology]
        float **topologyFreq;
        topologyFreq = new float *[nIter];
        if (topologyFreq == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyFreq[countI] = new float [nTopology];
	        if (topologyFreq == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq max [nIter]
        float *topologyFreqMax;
        topologyFreqMax = new float [nIter];
        if (topologyFreqMax == NULL)
            memoryError();

        // array to hold tree topology freq min [nIter]
        float *topologyFreqMin;
        topologyFreqMin = new float [nIter];
        if (topologyFreqMin == NULL)
            memoryError();

        // array to hold tree topology freq variance [nIter]
        float *topologyFreqVar;
        topologyFreqVar = new float [nIter];
        if (topologyFreqVar == NULL)
            memoryError();

    cout << "data arrays initiated" << endl;     
       
    // read input file
    readInputDsp(inputFileDsp, nDispStep, dispStep, dispProb);

	// dispInputSettingDsp
    dispInputSettingDsp(inputNameDsp, nPopSize, nDispStep, dispStep, dispProb);
    
    // getDispTable
    getDispTable(nDispStep, dispStep, dispProb, dispTable);
    
    // main loop structure
    // 1. iter loop (for)
    //   2. pedigree loop (for)
    //     3. generation loop (while)
    //       4. tree loop (for)   
    //         5. allele loop (for)

    // for each iter
    for (countI = 0; countI < nIter; countI++) {
        // reset event1, event2, tree topology
		for (countP = 0; countP < nPedigree; countP++) {
            for (countT = 0; countT < nTree; countT++) {
                event1[countP][countT] = 0;
                event2[countP][countT] = 0;
                treeTopology[countP][countT] = -1;
            }
		}
    
        // for each pedigree
        for (countP = 0; countP < nPedigree; countP++) {
            // reset preLoc, alleleCoal
            for (countT = 0; countT < nTree; countT++) {
                for (countA = 0; countA < nAllele; countA++) {
                    preLoc[countT][countA] = countA * distAllele;
                    alleleCoal[countT][countA] = countA; // coal with itself
                }
            }
            
            // reset generation count
            countG = 0;

            // reset flagRun
            flagRun = 1;

            // generation loop (while flagRun == true)
            while (flagRun == 1) {
                // update generation count
                countG++;

                // get p0Loc, p1Loc
                // for each individual
                for (countN = 0; countN < nPopSize; countN++) {
                    // countN is the location for the individual itself
                    p0Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                    p1Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                }  // for each individual

                // for each tree
                for (countT = 0; countT < nTree; countT++) {
                    // if event2 haven't occurred yet
                    if (event2[countP][countT] == 0) {
                        // for each allele
                        for (countA = 0; countA < nAllele; countA++) {
                            // update curLoc (copy from preLoc)
                            curLoc[countT][countA] = preLoc[countT][countA];

                            // get preLoc
                            // if coaled with itself (not coaled with another allele)
                            //   => get preLoc (randomly pick from p0 or p1)
                            // else (already coaled with another allele)
                            //   => no need to process (removed from pop), set preLoc = -1
                            if (alleleCoal[countT][countA] == countA){
                                // either p0 or p1
                                if (getRandomInt(0,1) == 0) {
                                    // if from p0
                                    preLoc[countT][countA] = p0Loc[(curLoc[countT][countA])];
                                }
                                else {
                                    // from p1
                                    preLoc[countT][countA] = p1Loc[(curLoc[countT][countA])];
                                }
                            }
                            else {
                                preLoc[countT][countA] = -1;
                            } // get preLoc
                        }   // for each allele
                        
                        // determine coal event
                        // for each allele pair
                        for (countA1 = 0; countA1 < (nAllele - 1); countA1++)  {
                            // only check with another allele if the allele exists
                            if (preLoc[countT][countA1] != -1) {
                                for (countA2 = (countA1 + 1); countA2 < nAllele; countA2++) {
                                    // if preLoc are the same => 1/2 prob of coal, else do nothing
                                    if (preLoc[countT][countA1] == preLoc[countT][countA2]) {
                                        // if coal event, else do nothing
                                        if (getRandomInt(0,1) == 1) {
                                            // check if event1/event2 occured
                                            if (event1[countP][countT] == 0) {
                                                // event1 not yet occured
                                                event1[countP][countT] = (countA1 * 10 + countA2);
                                            }
                                            else if (event2[countP][countT] == 0) {
                                                // event2 not yet occured
                                                event2[countP][countT] = (countA1 * 10 + countA2);
                                            } // else = both event1 and event2 occurred, do nothing
                                            // The case for else here:
                                            // usually will not happen, 
                                            // but there's a rare chance that 
                                            // the 3 remaining alleles all coal in the same generation
                                            
                                            // A2 coaled with A1
                                            alleleCoal[countT][countA2] = countA1; 
                                            // remove A2
                                            preLoc[countT][countA2] = -1;
                                        } // if coal event, else do nothing
                                    } // if preLoc are the same, else do nothing
                                } // for countA2
                            } // only check with another allele if the allele exists, else do nothing
                        } // for countA1
                        // for each allele pair
                        
                    }// if event2 haven't occurred yet, else (event2 known) do nothing
                }   // for each tree
                
                // update flagRun
                flagRun = 0; // set to false (terminate generation loop)
                for (countT = 0; countT < nTree; countT++) {
                    // if find any undetermined event2
                    if (event2[countP][countT] == 0) {
                        flagRun = 1; // set flagRun to true
                        break;  // no need to check others
                    }
                } // update flagRun
            }   // generation loop (while flagRun == true)
        }    // for each pedigree    
        
        // getTreeTopology4
        getTreeTopology4(nPedigree, nTree, event1, event2, treeTopology);
        
        // get topologyCount, topologyFreq
        // statTopology
        statTopology(countI, nTopology, nPedigree, nTree, treeTopology, topologyCount, topologyFreq);
        
        // update topologyFreqMax, topologyFreqMin 
        // statTopologyFreq
        statTopologyFreq(countI, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

        // disp
        if (countI % repInterval == 0) {
            cout << endl;
            cout << "\t" << countI << "\t";
        }
        cout << '*';
    }    // for each iter
    
    // write to output
    // writeSummary
    writeSummary(nIter, topologyFreqMax, topologyFreqMin, topologyFreqVar);
    
    // writeTopologyFreq
    writeTopologyFreq(nIter, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

	// writeTopologyVar
	writeTopologyVar(outputFileRep, nIter, topologyFreqVar);

} // module1

// module2
// 4 alleles, unequal spacing at [0, 0.125, 0.25, 0.5]
void module2(void) {
    // nAllele
	const int nAllele = 4;
    // nTopology
	const int nTopology = 18;

    // file operation
        // input file - dispersal prob file
            // declare fileName
            fileName inputNameDsp; 
            // declare input file
            ifstream inputFileDsp;
            // call openInputFile function
            openInputFileDsp(inputFileDsp, inputNameDsp);

        // report file 
            // declare fileName
            fileName outputNameRep; 
            // declare output file
            ofstream outputFileRep;
            // call openOutput function
            openOutputRep(outputFileRep, outputNameRep);

        //
        writeHeader();

    // declare variables 
        // settings
        int nPopSize;   // pop size
        int nDispStep;  // number of possible dispersal steps 

        int nIter;      // number of iterations
        int nPedigree;  // number of pedigrees
        int nTree;      // number of gene trees per pedigree

        // variable
        int distAllele; // distance between alleles 
        int flagRun;    // bool flag for while loop

        // loop counting
        int countI;     // count iteration
        int countG;     // count generation
        int countN;     // count individual
        int countP;     // count pedigree
        int countT;     // count tree
        int countA;     // count allele
        int countA1;    // count allele
        int countA2;    // count allele

    // get settings from input file
    // readInputSettingDsp
    readInputSettingDsp(inputFileDsp, nPopSize, nDispStep);

    // distAllele (set to equal dist)
    distAllele = nPopSize / nAllele;
            
    // get settings from user
    getNIter(nIter);
    getNPedigree(nPedigree);
    getNTree(nTree);
    
    // write settings to stdout
    writeSetting(inputNameDsp, outputNameRep, nPopSize, nDispStep, nIter, nPedigree, nTree);
    
    // declare and initiate array   
    cout << "initiating data arrays..." << endl;     
        // array to hold dispStep
        int *dispStep;
        dispStep = new int [nDispStep];
        if (dispStep == NULL)
            memoryError();        
        
        // array to hold dispProb 
        float *dispProb;
        dispProb = new float [nDispStep];
        if (dispProb == NULL)
            memoryError();        

        // array to hold dispTable
        // hush table to find nStep
        int *dispTable;
        dispTable = new int [(maxRandomInt + 1)];
        if (dispTable == NULL)
            memoryError();
        
        // parent location
        int *p0Loc;
        p0Loc = new int [nPopSize];
        if (p0Loc == NULL)
            memoryError();
        
        // parent location
        int *p1Loc;
        p1Loc = new int [nPopSize];
        if (p1Loc == NULL)
            memoryError();
                    
        // array to hold location in current generation [nTree][nAllele]
        int **curLoc;
        curLoc = new int *[nTree];
        if (curLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			curLoc[countT] = new int [nAllele];
	        if (curLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold location in previous generation [nTree][nAllele]
        int **preLoc;
        preLoc = new int *[nTree];
        if (preLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			preLoc[countT] = new int [nAllele];
	        if (preLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold allele that has been coaled [nTree][nAllele]
        // either itself (have not coaled with any other allele)
        // or a smaller alleleId
        // when more than 2 alleles coaled, use the smallest
        //   e.g if allele 1, 2, 4 coaled
        //   allele[countT][0] = 0
        //   allele[countT][1] = 0
        //   allele[countT][3] = 0
        int **alleleCoal;
        alleleCoal = new int *[nTree];
        if (alleleCoal == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			alleleCoal[countT] = new int [nAllele];
	        if (alleleCoal == NULL)
	            memoryError();			 
		}
        
        // array to hold coal event1 id [nPedigree][nTree]
        int **event1;
        event1 = new int *[nPedigree];
        if (event1 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event1[countP] = new int [nTree];
	        if (event1 == NULL)
	            memoryError();			 
		}

        // array to hold coal event2 id [nPedigree][nTree]
        int **event2;
        event2 = new int *[nPedigree];
        if (event2 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event2[countP] = new int [nTree];
	        if (event2 == NULL)
	            memoryError();			 
		}

        // array to hold tree topology id [nPedigree][nTree]
        int **treeTopology;
        treeTopology = new int *[nPedigree];
        if (treeTopology == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			treeTopology[countP] = new int [nTree];
	        if (treeTopology == NULL)
	            memoryError();			 
		}
        
        // data summary arrays
        
        // array to hold tree topology count [nIter][nTopology]
        int **topologyCount;
        topologyCount = new int *[nIter];
        if (topologyCount == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyCount[countI] = new int [nTopology];
	        if (topologyCount == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq [nIter][nTopology]
        float **topologyFreq;
        topologyFreq = new float *[nIter];
        if (topologyFreq == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyFreq[countI] = new float [nTopology];
	        if (topologyFreq == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq max [nIter]
        float *topologyFreqMax;
        topologyFreqMax = new float [nIter];
        if (topologyFreqMax == NULL)
            memoryError();

        // array to hold tree topology freq min [nIter]
        float *topologyFreqMin;
        topologyFreqMin = new float [nIter];
        if (topologyFreqMin == NULL)
            memoryError();

        // array to hold tree topology freq variance [nIter]
        float *topologyFreqVar;
        topologyFreqVar = new float [nIter];
        if (topologyFreqVar == NULL)
            memoryError();

    cout << "data arrays initiated" << endl;     
       
    // read input file
    readInputDsp(inputFileDsp, nDispStep, dispStep, dispProb);

	// dispInputSettingDsp
    dispInputSettingDsp(inputNameDsp, nPopSize, nDispStep, dispStep, dispProb);
    
    // getDispTable
    getDispTable(nDispStep, dispStep, dispProb, dispTable);
    
    // main loop structure
    // 1. iter loop (for)
    //   2. pedigree loop (for)
    //     3. generation loop (while)
    //       4. tree loop (for)   
    //         5. allele loop (for)

    // for each iter
    for (countI = 0; countI < nIter; countI++) {
        // reset event1, event2, tree topology
		for (countP = 0; countP < nPedigree; countP++) {
            for (countT = 0; countT < nTree; countT++) {
                event1[countP][countT] = 0;
                event2[countP][countT] = 0;
                treeTopology[countP][countT] = -1;
            }
		}
    
        // for each pedigree
        for (countP = 0; countP < nPedigree; countP++) {
            // reset preLoc, alleleCoal
            for (countT = 0; countT < nTree; countT++) {
            	// 4 alleles, unequal spacing at [0, 0.125, 0.25, 0.5]
                preLoc[countT][0] = 0;
                preLoc[countT][1] = int (0.125 * nPopSize);
                preLoc[countT][2] = int (0.25 * nPopSize);
                preLoc[countT][3] = int (0.5 * nPopSize);
                for (countA = 0; countA < nAllele; countA++) {
                    alleleCoal[countT][countA] = countA; // coal with itself
                }
            }
            
            // reset generation count
            countG = 0;

            // reset flagRun
            flagRun = 1;

            // generation loop (while flagRun == true)
            while (flagRun == 1) {
                // update generation count
                countG++;

                // get p0Loc, p1Loc
                // for each individual
                for (countN = 0; countN < nPopSize; countN++) {
                    // countN is the location for the individual itself
                    p0Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                    p1Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                }  // for each individual

                // for each tree
                for (countT = 0; countT < nTree; countT++) {
                    // if event2 haven't occurred yet
                    if (event2[countP][countT] == 0) {
                        // for each allele
                        for (countA = 0; countA < nAllele; countA++) {
                            // update curLoc (copy from preLoc)
                            curLoc[countT][countA] = preLoc[countT][countA];

                            // get preLoc
                            // if coaled with itself (not coaled with another allele)
                            //   => get preLoc (randomly pick from p0 or p1)
                            // else (already coaled with another allele)
                            //   => no need to process (removed from pop), set preLoc = -1
                            if (alleleCoal[countT][countA] == countA){
                                // either p0 or p1
                                if (getRandomInt(0,1) == 0) {
                                    // if from p0
                                    preLoc[countT][countA] = p0Loc[(curLoc[countT][countA])];
                                }
                                else {
                                    // from p1
                                    preLoc[countT][countA] = p1Loc[(curLoc[countT][countA])];
                                }
                            }
                            else {
                                preLoc[countT][countA] = -1;
                            } // get preLoc
                        }   // for each allele
                        
                        // determine coal event
                        // for each allele pair
                        for (countA1 = 0; countA1 < (nAllele - 1); countA1++)  {
                            // only check with another allele if the allele exists
                            if (preLoc[countT][countA1] != -1) {
                                for (countA2 = (countA1 + 1); countA2 < nAllele; countA2++) {
                                    // if preLoc are the same => 1/2 prob of coal, else do nothing
                                    if (preLoc[countT][countA1] == preLoc[countT][countA2]) {
                                        // if coal event, else do nothing
                                        if (getRandomInt(0,1) == 1) {
                                            // check if event1/event2 occured
                                            if (event1[countP][countT] == 0) {
                                                // event1 not yet occured
                                                event1[countP][countT] = (countA1 * 10 + countA2);
                                            }
                                            else if (event2[countP][countT] == 0) {
                                                // event2 not yet occured
                                                event2[countP][countT] = (countA1 * 10 + countA2);
                                            } // else = both event1 and event2 occurred, do nothing
                                            // The case for else here:
                                            // usually will not happen, 
                                            // but there's a rare chance that 
                                            // the 3 remaining alleles all coal in the same generation
                                            
                                            // A2 coaled with A1
                                            alleleCoal[countT][countA2] = countA1; 
                                            // remove A2
                                            preLoc[countT][countA2] = -1;
                                        } // if coal event, else do nothing
                                    } // if preLoc are the same, else do nothing
                                } // for countA2
                            } // only check with another allele if the allele exists, else do nothing
                        } // for countA1
                        // for each allele pair
                        
                    }// if event2 haven't occurred yet, else (event2 known) do nothing
                }   // for each tree
                
                // update flagRun
                flagRun = 0; // set to false (terminate generation loop)
                for (countT = 0; countT < nTree; countT++) {
                    // if find any undetermined event2
                    if (event2[countP][countT] == 0) {
                        flagRun = 1; // set flagRun to true
                        break;  // no need to check others
                    }
                } // update flagRun
            }   // generation loop (while flagRun == true)
        }    // for each pedigree    
        
        // getTreeTopology4
        getTreeTopology4(nPedigree, nTree, event1, event2, treeTopology);
        
        // get topologyCount, topologyFreq
        // statTopology
        statTopology(countI, nTopology, nPedigree, nTree, treeTopology, topologyCount, topologyFreq);
        
        // update topologyFreqMax, topologyFreqMin 
        // statTopologyFreq
        statTopologyFreq(countI, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

        // disp
        if (countI % repInterval == 0) {
            cout << endl;
            cout << "\t" << countI << "\t";
        }
        cout << '*';
    }    // for each iter
    
    // write to output
    // writeSummary
    writeSummary(nIter, topologyFreqMax, topologyFreqMin, topologyFreqVar);
    
    // writeTopologyFreq
    writeTopologyFreq(nIter, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

	// writeTopologyVar
	writeTopologyVar(outputFileRep, nIter, topologyFreqVar);

} // module2

// module3
// 5 alleles, equal spacing at [0, 0.2, 0.4, 0.6, 0.8]
void module3(void) {
    // nAllele
	const int nAllele = 5;
    // nTopology
	const int nTopology = 180;

    // file operation
        // input file - dispersal prob file
            // declare fileName
            fileName inputNameDsp; 
            // declare input file
            ifstream inputFileDsp;
            // call openInputFile function
            openInputFileDsp(inputFileDsp, inputNameDsp);

        // report file 
            // declare fileName
            fileName outputNameRep; 
            // declare output file
            ofstream outputFileRep;
            // call openOutput function
            openOutputRep(outputFileRep, outputNameRep);

        //
        writeHeader();

    // declare variables 
        // settings
        int nPopSize;   // pop size
        int nDispStep;  // number of possible dispersal steps 

        int nIter;      // number of iterations
        int nPedigree;  // number of pedigrees
        int nTree;      // number of gene trees per pedigree

        // variable
        int distAllele; // distance between alleles (= nPopSize / 4)
        int flagRun;    // bool flag for while loop

        // loop counting
        int countI;     // count iteration
        int countG;     // count generation
        int countN;     // count individual
        int countP;     // count pedigree
        int countT;     // count tree
        int countA;     // count allele
        int countA1;    // count allele
        int countA2;    // count allele

    // get settings from input file
    // readInputSettingDsp
    readInputSettingDsp(inputFileDsp, nPopSize, nDispStep);

    // distAllele (set to equal dist)
    distAllele = nPopSize / nAllele;
            
    // get settings from user
    getNIter(nIter);
    getNPedigree(nPedigree);
    getNTree(nTree);
    
    // write settings to stdout
    writeSetting(inputNameDsp, outputNameRep, nPopSize, nDispStep, nIter, nPedigree, nTree);
    
    // declare and initiate array   
    cout << "initiating data arrays..." << endl;     
        // array to hold dispStep
        int *dispStep;
        dispStep = new int [nDispStep];
        if (dispStep == NULL)
            memoryError();        
        
        // array to hold dispProb 
        float *dispProb;
        dispProb = new float [nDispStep];
        if (dispProb == NULL)
            memoryError();        

        // array to hold dispTable
        // hush table to find nStep
        int *dispTable;
        dispTable = new int [(maxRandomInt + 1)];
        if (dispTable == NULL)
            memoryError();
        
        // parent location
        int *p0Loc;
        p0Loc = new int [nPopSize];
        if (p0Loc == NULL)
            memoryError();
        
        // parent location
        int *p1Loc;
        p1Loc = new int [nPopSize];
        if (p1Loc == NULL)
            memoryError();
                    
        // array to hold location in current generation [nTree][nAllele]
        int **curLoc;
        curLoc = new int *[nTree];
        if (curLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			curLoc[countT] = new int [nAllele];
	        if (curLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold location in previous generation [nTree][nAllele]
        int **preLoc;
        preLoc = new int *[nTree];
        if (preLoc == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			preLoc[countT] = new int [nAllele];
	        if (preLoc == NULL)
	            memoryError();			 
		}
        
        // array to hold allele that has been coaled [nTree][nAllele]
        // either itself (have not coaled with any other allele)
        // or a smaller alleleId
        // when more than 2 alleles coaled, use the smallest
        //   e.g if allele 1, 2, 4 coaled
        //   allele[countT][0] = 0
        //   allele[countT][1] = 0
        //   allele[countT][3] = 0
        int **alleleCoal;
        alleleCoal = new int *[nTree];
        if (alleleCoal == NULL)
            memoryError();

		for (countT = 0; countT < nTree; countT++) {
			alleleCoal[countT] = new int [nAllele];
	        if (alleleCoal == NULL)
	            memoryError();			 
		}
        
        // array to hold coal event1 id [nPedigree][nTree]
        int **event1;
        event1 = new int *[nPedigree];
        if (event1 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event1[countP] = new int [nTree];
	        if (event1 == NULL)
	            memoryError();			 
		}

        // array to hold coal event2 id [nPedigree][nTree]
        int **event2;
        event2 = new int *[nPedigree];
        if (event2 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event2[countP] = new int [nTree];
	        if (event2 == NULL)
	            memoryError();			 
		}

        // array to hold coal event3 id [nPedigree][nTree]
        int **event3;
        event3 = new int *[nPedigree];
        if (event3 == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			event3[countP] = new int [nTree];
	        if (event3 == NULL)
	            memoryError();			 
		}

        // array to hold tree topology id [nPedigree][nTree]
        int **treeTopology;
        treeTopology = new int *[nPedigree];
        if (treeTopology == NULL)
            memoryError();

		for (countP = 0; countP < nPedigree; countP++) {
			treeTopology[countP] = new int [nTree];
	        if (treeTopology == NULL)
	            memoryError();			 
		}
        
        // data summary arrays
        
        // array to hold tree topology count [nIter][nTopology]
        int **topologyCount;
        topologyCount = new int *[nIter];
        if (topologyCount == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyCount[countI] = new int [nTopology];
	        if (topologyCount == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq [nIter][nTopology]
        float **topologyFreq;
        topologyFreq = new float *[nIter];
        if (topologyFreq == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++) {
			topologyFreq[countI] = new float [nTopology];
	        if (topologyFreq == NULL)
	            memoryError();			 
		}
        
        // array to hold tree topology freq max [nIter]
        float *topologyFreqMax;
        topologyFreqMax = new float [nIter];
        if (topologyFreqMax == NULL)
            memoryError();

        // array to hold tree topology freq min [nIter]
        float *topologyFreqMin;
        topologyFreqMin = new float [nIter];
        if (topologyFreqMin == NULL)
            memoryError();

        // array to hold tree topology freq variance [nIter]
        float *topologyFreqVar;
        topologyFreqVar = new float [nIter];
        if (topologyFreqVar == NULL)
            memoryError();

    cout << "data arrays initiated" << endl;     
       
    // read input file
    readInputDsp(inputFileDsp, nDispStep, dispStep, dispProb);

	// dispInputSettingDsp
    dispInputSettingDsp(inputNameDsp, nPopSize, nDispStep, dispStep, dispProb);
    
    // getDispTable
    getDispTable(nDispStep, dispStep, dispProb, dispTable);
    
    // main loop structure
    // 1. iter loop (for)
    //   2. pedigree loop (for)
    //     3. generation loop (while)
    //       4. tree loop (for)   
    //         5. allele loop (for)

    // for each iter
    for (countI = 0; countI < nIter; countI++) {
        // reset event1, event2, event3, tree topology
		for (countP = 0; countP < nPedigree; countP++) {
            for (countT = 0; countT < nTree; countT++) {
                event1[countP][countT] = 0;
                event2[countP][countT] = 0;
                event3[countP][countT] = 0;
                treeTopology[countP][countT] = -1;
            }
		}
    
        // for each pedigree
        for (countP = 0; countP < nPedigree; countP++) {
            // reset preLoc, alleleCoal
            for (countT = 0; countT < nTree; countT++) {
                for (countA = 0; countA < nAllele; countA++) {
                    preLoc[countT][countA] = countA * distAllele;
                    alleleCoal[countT][countA] = countA; // coal with itself
                }
            }
            
            // reset generation count
            countG = 0;

            // reset flagRun
            flagRun = 1;

            // generation loop (while flagRun == true)
            while (flagRun == 1) {
                // update generation count
                countG++;

                // get p0Loc, p1Loc
                // for each individual
                for (countN = 0; countN < nPopSize; countN++) {
                    // countN is the location for the individual itself
                    p0Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                    p1Loc[countN] = getParLoc(countN, nPopSize, dispTable);
                }  // for each individual

                // for each tree
                for (countT = 0; countT < nTree; countT++) {
                    // if event3 haven't occurred yet
                    if (event3[countP][countT] == 0) {
                        // for each allele
                        for (countA = 0; countA < nAllele; countA++) {
                            // update curLoc (copy from preLoc)
                            curLoc[countT][countA] = preLoc[countT][countA];

                            // get preLoc
                            // if coaled with itself (not coaled with another allele)
                            //   => get preLoc (randomly pick from p0 or p1)
                            // else (already coaled with another allele)
                            //   => no need to process (removed from pop), set preLoc = -1
                            if (alleleCoal[countT][countA] == countA){
                                // either p0 or p1
                                if (getRandomInt(0,1) == 0) {
                                    // if from p0
                                    preLoc[countT][countA] = p0Loc[(curLoc[countT][countA])];
                                }
                                else {
                                    // from p1
                                    preLoc[countT][countA] = p1Loc[(curLoc[countT][countA])];
                                }
                            }
                            else {
                                preLoc[countT][countA] = -1;
                            } // get preLoc
                        }   // for each allele
                        
                        // determine coal event
                        // for each allele pair
                        for (countA1 = 0; countA1 < (nAllele - 1); countA1++)  {
                            // only check with another allele if the allele exists
                            if (preLoc[countT][countA1] != -1) {
                                for (countA2 = (countA1 + 1); countA2 < nAllele; countA2++) {
                                    // if preLoc are the same => 1/2 prob of coal, else do nothing
                                    if (preLoc[countT][countA1] == preLoc[countT][countA2]) {
                                        // if coal event, else do nothing
                                        if (getRandomInt(0,1) == 1) {
                                            // check if event1/event2 occured
                                            if (event1[countP][countT] == 0) {
                                                // event1 not yet occured
                                                event1[countP][countT] = (countA1 * 10 + countA2);
                                            }
                                            else if (event2[countP][countT] == 0) {
                                                // event2 not yet occured
                                                event2[countP][countT] = (countA1 * 10 + countA2);
                                            } 
                                            else if (event3[countP][countT] == 0) {
                                                // event3 not yet occured
                                                event3[countP][countT] = (countA1 * 10 + countA2);
                                            } 
                                            else {
                                            	// else = all 3 events occurred, do nothing
                                            }
                                            // A2 coaled with A1
                                            alleleCoal[countT][countA2] = countA1; 
                                            // remove A2
                                            preLoc[countT][countA2] = -1;
                                        } // if coal event, else do nothing
                                    } // if preLoc are the same, else do nothing
                                } // for countA2
                            } // only check with another allele if the allele exists, else do nothing
                        } // for countA1
                        // for each allele pair
                        
                    }// if event3 haven't occurred yet, else (event3 known) do nothing
                }   // for each tree
                
                // update flagRun
                flagRun = 0; // set to false (terminate generation loop)
                for (countT = 0; countT < nTree; countT++) {
                    // if find any undetermined event3
                    if (event3[countP][countT] == 0) {
                        flagRun = 1; // set flagRun to true
                        break;  // no need to check others
                    }
                } // update flagRun
            }   // generation loop (while flagRun == true)
        }    // for each pedigree    
        
        // getTreeTopology5
        getTreeTopology5(nPedigree, nTree, event1, event2, event3, treeTopology);
        
        // get topologyCount, topologyFreq
        // statTopology
        statTopology(countI, nTopology, nPedigree, nTree, treeTopology, topologyCount, topologyFreq);
        
        // update topologyFreqMax, topologyFreqMin 
        // statTopologyFreq
        statTopologyFreq(countI, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

        // disp
        if (countI % repInterval == 0) {
            cout << endl;
            cout << "\t" << countI << "\t";
        }
        cout << '*';
    }    // for each iter
    
    // write to output
    // writeSummary
    writeSummary(nIter, topologyFreqMax, topologyFreqMin, topologyFreqVar);
    
    // writeTopologyFreq
    writeTopologyFreq(nIter, nTopology, topologyFreq, topologyFreqMax, topologyFreqMin, topologyFreqVar);

	// writeTopologyVar
	writeTopologyVar(outputFileRep, nIter, topologyFreqVar);

} // module3

// main function
int main(void) {
    // call welcome function to display welcome message
    welcome();
        
    // seed random number generator with current system time
    srand(time(0));

    // declare variables
    int module; //module used for the run

    do {        
        getModule(module);
        // call selected module
        switch (module) {
            case 0: break;
            case 1: module1();
                    break;
            case 2: module2();
                    break;
            case 3: module3();
                    break;
            default: cout << endl;
                     cout << "*** Module selection error ***" << endl;  
                     exit(0);
        } // end of switch    
    } while (module != 0);

	// end, display message
    	cout << endl;
    	cout << "*** Program ended. ***" << endl;
        return 0;   
} // end of main function