// phyloBreak source code
// Chih-Horng Kuo
// chkuo@lifedev.org


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

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

// constants 
const char programName[21] = "phyloBreak";
    // program name
const char versionNumber[15] = "A2";
    // version number
const char versionDate[15] = "04/12/2004";
    // version date
const char author[61] = "Chih-Horng Kuo";
    // author
const char email[21] = "chkuo@lifedev.org";
    // author email
const char web[61] = "http://chkuo.name/";
    // author web
const int fileNameLength = 81;
    // filename length
const int maxRandomInt = 20000;
    // max random int for generating random float    
const char symbolComment = '/';
    // symbol in input/output file for comment    
const char symbolNIter = '#';
    // symbol in input/output file for nIter    
const char symbolOptSize = '*';
    // symbol in input/output file for optSize    
const char symbolOptBoundary = '|';
    // symbol in input/output file for optBoundary    
const char symbolNDispStep = '^';
    // symbol in input/output file for nDispStep    
const char symbolDispProb = '@';
    // symbol in input/output file for the start of dispersal prob table
const char symbolMCT = '$';
    // symbol in input/output file for the start of coalescent time matrix
const char symbolMCA = '!';
    // symbol in input/output file for the start of coefficient of association matrix
const char symbolMGD = '%';
    // symbol in input/output file for the start of average genetic distance matrix
    

// typedef
typedef char fileName[fileNameLength];
    // new data type "fileName" to hold file name

// function prototypes
    // infomation display
    void welcome(void);
    void memoryError(void);

    // get option
    void getModule(int &module);

    // modules
    void module1(void);


// start of program functions
// 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 1: module1();
                        break;
                case 2: break;
                default: cout << endl;
                         cout << "*** Module selection error" << endl;  
                         exit(0);
            } // end of switch    
    } while (module != 2);

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



// welcome function
// display welcome info
void welcome(void)
{
    // display information for users
    cout << "*** " << programName << " ***" << endl << endl;
    cout << "Version number: " << versionNumber << endl;
    cout << "Release date: " << versionDate << endl;
    cout << endl;

    cout << "Contact information" << endl; 
    cout << "  Author: " << author << endl; 
    cout << "  email: " << email << endl;
    cout << "  Web: " << web << endl;
    cout << endl;

    cout << "Description:" << 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


// openOutputLog function
// open output file and get output file name
void openOutputLog(ofstream &outputFileLog, fileName &outputNameLog)
{
    // 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 (phyloBreak Log <.pbl>): ";
    cin >> outputNameLog;
    
    // open outputFile
        outputFileLog.open(outputNameLog);
        if ( outputFileLog == 0 )
        {
    	    cout << endl;
    		cout << "*** Error opening output file." << endl;
            exit(0);
        }
        else
		cout << endl;
		cout << "Output file (" << outputNameLog << ") opened successfully." << endl;

} // end of openOutputLog


// openOutputMct function
// open output file and get output file name
void openOutputMct(ofstream &outputFileMct, fileName &outputNameMct)
{
    // 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 (Coalescent time matrix <.mct>): ";
    cin >> outputNameMct;
    
    // open outputFile
        outputFileMct.open(outputNameMct);
        if ( outputFileMct == 0 )
        {
    	    cout << endl;
    		cout << "*** Error opening output file." << endl;
            exit(0);
        }
        else
		cout << endl;
		cout << "Output file (" << outputNameMct << ") opened successfully." << endl;

} // end of openOutputMct

// openOutputMca function
// open output file and get output file name
void openOutputMca(ofstream &outputFileMca, fileName &outputNameMca)
{
    // 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 (Coefficient of association matrix <.mca>): ";
    cin >> outputNameMca;
    
    // open outputFile
        outputFileMca.open(outputNameMca);
        if ( outputFileMca == 0 )
        {
    	    cout << endl;
    		cout << "*** Error opening output file." << endl;
            exit(0);
        }
        else
		cout << endl;
		cout << "Output file (" << outputNameMca << ") opened successfully." << endl;

} // end of openOutputMca

// openOutputMgd function
// open output file and get output file name
void openOutputMgd(ofstream &outputFileMgd, fileName &outputNameMgd)
{
    // 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 (Average genetic distance matrix <.mgd>): ";
    cin >> outputNameMgd;
    
    // open outputFile
        outputFileMgd.open(outputNameMgd);
        if ( outputFileMgd == 0 )
        {
    	    cout << endl;
    		cout << "*** Error opening output file." << endl;
            exit(0);
        }
        else
		cout << endl;
		cout << "Output file (" << outputNameMgd << ") opened successfully." << endl;

} // end of openOutputMgd



// closeInputFile function
// close input file
void closeInputFile(ifstream &inputFile)
{
    cout << endl;
    
    inputFile.close();
}


// closeOutputFile function
// close output file
void closeOutputFile(ofstream &outputFile)
{
    cout << endl;
    
    outputFile.close();
}


// readInputSettingDsp function
void readInputSettingDsp(ifstream &inputFileDsp, int &nIter, int &nPopSize, int &nLoc, int &nInd, int &optBoundary, int &posBarrier, int &durBarrier, int &nDispStep)
{
    // declare variables
    char symbolInput;
    char inputComment[81];
    int flagNIter = 0;
    int flagOptSize = 0;
    int flagOptBoundary = 0;
    int flagNDispStep = 0;
    
    // read input setting
    do
    {
        inputFileDsp.get(symbolInput);
        switch (symbolInput)
        {
            case symbolComment:     inputFileDsp.getline(inputComment, 81, '\n');
                                    break;
            case symbolNIter:       inputFileDsp >> nIter;
                                    flagNIter = 1;
                                    break;
            case symbolOptSize:     inputFileDsp >> nPopSize;
                                    inputFileDsp >> nLoc;
                                    inputFileDsp >> nInd;                                    
                                    flagOptSize = 1;
                                    break;
            case symbolOptBoundary: inputFileDsp >> optBoundary;
                                    inputFileDsp >> posBarrier;
                                    inputFileDsp >> durBarrier;
                                    flagOptBoundary = 1;
                                    break;
            case symbolNDispStep:   inputFileDsp >> nDispStep;
                                    flagNDispStep = 1;
                                    break;
        }
    } while (flagNIter !=1 || flagOptSize !=1 || flagOptBoundary !=1 || flagNDispStep !=1);

    // check input setting

    // nIter
	if (nIter < 1)
	{
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** nIter < 1" << endl;
		exit(0);
	}    
    // optSize
    // check if nPopSize < nLoc*nInd
    if (nPopSize < nLoc*nInd)
    {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** nPopSize < nLoc*nInd" << endl;
		exit(0);
    }


    // optBoundary
    // 1 = No boundaries (donut-shape habitat)
    // 2 = Absorbing boundaries
    // 3 = Reflecting boundaries
    if (optBoundary < 1 || optBoundary > 3)
    {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** optBoundary setting error" << endl;
		exit(0);
    }
    
    if (posBarrier < 0 || posBarrier > nPopSize)
    {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** posBarrier setting error" << endl;
		exit(0);
    }
    
    if (durBarrier < 0)
    {
	    cout << endl;
		cout << "*** Input setting error" << endl;
		cout << "*** durBarrier setting error" << 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
    //
} // readInputDsp

// dispInputSettingDsp
void dispInputSettingDsp(fileName &inputNameDsp, int nIter, int nPopSize, int nLoc, int nInd, int optBoundary, int posBarrier, int durBarrier, int nDispStep, int *dispStep, float *dispProb)
{
	//
	int countD;
	
	//
	cout << endl;
	cout << "Input setting from " << inputNameDsp << endl;
	cout << endl;
	
	cout << "  nIter = " << nIter << endl;
	cout << "  nPopSize = " << nPopSize << endl;
	cout << "  nLoc = " << nLoc << endl;
	cout << "  nInd = " << nInd << endl;
	cout << "  optBoundary = " << optBoundary << endl;
	cout << "  posBarrier = " << posBarrier << endl;
	cout << "  durBarrier = " << durBarrier << endl;
	cout << "  nDispStep = " << nDispStep << endl;
	cout << endl;
	
	cout << "  dispStep" << "\t" << "dispProb" << endl;
	for (countD = 0; countD < nDispStep; countD++)
	{
		cout << "  " << dispStep[countD] << "\t" << dispProb[countD] << endl;
	}
	cout << endl;	

} // dispInputSettingDsp


// writeOutputHeaderM1 function
void writeOutputHeaderM1(ofstream &outputFileLog, fileName outputNameLog)
{   
    cout << "Writing header to output files..." << endl;
    cout << endl;

    // get time
    time_t rawtime;
    struct tm * timeinfo;

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

    // write header to Log file
    outputFileLog << " *** Start of header ***" << endl;
    outputFileLog << " Program name: " << programName << endl;
    outputFileLog << "   Version number: " << versionNumber << endl;
    outputFileLog << "   Release date: " << versionDate << endl;
    outputFileLog << endl;

    outputFileLog << " Contact information" << endl; 
    outputFileLog << "   Author: " << author << endl; 
    outputFileLog << "   email: " << email << endl;
    outputFileLog << "   Web: " << web << endl;
    outputFileLog << endl;
    
    outputFileLog << " Run date and time: " << asctime(timeinfo);
    outputFileLog << " Filename: " << outputNameLog << endl;
    outputFileLog << " *** End of header ***" << endl;
    outputFileLog << endl;

} // end of writeOutputHeaderM1 function

// writeM1Setting function
void writeM1Setting(ofstream &outputFileLog, fileName outputNameLog, fileName outputNameMct, fileName outputNameMca, fileName outputNameMgd, fileName inputNameDsp, int nIter, int nPopSize, int nLoc, int nInd, int optBoundary, int posBarrier, int durBarrier, int nDispStep)
{
    // write settings to Log
    outputFileLog << " *** Start of settings ***" << endl;
    outputFileLog << " Input file" << endl;
    outputFileLog << "   Dsp file: " << inputNameDsp << endl;
    outputFileLog << " Output file" << endl;
    outputFileLog << "   Log file: " << outputNameLog << endl;
    outputFileLog << "   Mct file: " << outputNameMct << endl;
    outputFileLog << "   Mca file: " << outputNameMca << endl;
    outputFileLog << "   Mgd file: " << outputNameMgd << endl;
    outputFileLog << endl;

    outputFileLog << " input settings" << endl; 
    outputFileLog << "   nIter: " << nIter << endl;
    outputFileLog << "   nPopSize: " << nPopSize << endl;
    outputFileLog << "   nLoc: " << nLoc << endl;
    outputFileLog << "   nInd: " << nInd << endl;
    outputFileLog << "   optBoundary: " << optBoundary << endl;
    outputFileLog << "   posBarrier: " << posBarrier << endl;
    outputFileLog << "   durBarrier: " << durBarrier << endl;
    outputFileLog << "   nDispStep: " << nDispStep << endl;
    outputFileLog << endl;
    
    outputFileLog << " *** End of settings ***" << endl;
    outputFileLog << endl;
} // writeM1Setting function

// writeOutputSetting
void writeOutputSetting(ofstream &outputFile, int nIter, int nPopSize, int nLoc, int nInd, int optBoundary, int posBarrier, int durBarrier, int nDispStep)
{

	outputFile << symbolNIter << "\t" << nIter << endl;
	outputFile << endl;

	outputFile << symbolOptSize << "\t" 
			   << nPopSize << "\t"
			   << nLoc << "\t"
			   << nInd << endl;
	outputFile << endl;

	outputFile << symbolOptBoundary << "\t" 
			   << optBoundary << "\t"
			   << posBarrier << "\t"
			   << durBarrier << endl;
	outputFile << endl;
	
	outputFile << symbolNDispStep << "\t" << nDispStep << endl;
	outputFile << endl;

} // writeOutputSetting


// writeMct
void writeMct(ofstream &outputFileMct, int countI, int countG, int nSample, int **coalTime)
{
    // 
    int countRow;
    int countCol;
    
    outputFileMct << symbolMCT << "\t"
                  << countI << "\t"
                  << countG << endl;
    
    for (countRow = 0; countRow < nSample; countRow++)
    {
        for (countCol = 0; countCol < nSample; countCol++)
        {
            outputFileMct << coalTime[countRow][countCol] << "\t";
        }
        outputFileMct << endl;
    }
    

} // writeMct

// writeMca
void writeMca(ofstream &outputFileMca, int countI, int countG, int nLoc, float **arrayCOA)
{
    // 
    int countRow;
    int countCol;
    
    outputFileMca << symbolMCA << "\t"
                  << countI << "\t"
                  << countG << endl;
    
    for (countRow = 0; countRow < nLoc; countRow++)
    {
        for (countCol = 0; countCol < nLoc; countCol++)
        {
            outputFileMca << setprecision(4)
	                      << setiosflags(ios::fixed | ios::showpoint)
    	                  << arrayCOA[countRow][countCol] << "\t";
        }
        outputFileMca << endl;
    }    

} // writeMca

// writeMgd
void writeMgd(ofstream &outputFileMgd, int countI, int countG, int nLoc, float **avGD)
{
    // 
    int countRow;
    int countCol;
    
    outputFileMgd << symbolMGD << "\t"
                  << countI << "\t"
                  << countG << endl;
    
    for (countRow = 0; countRow < nLoc; countRow++)
    {
        for (countCol = 0; countCol < nLoc; countCol++)
        {
            outputFileMgd << setprecision(4)
	                      << setiosflags(ios::fixed | ios::showpoint)
    	                  << avGD[countRow][countCol] << "\t";
        }
        outputFileMgd << endl;
    }    

} // writeMgd


// 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 << "Pedigree generator modules available:" << endl;
    cout << "  (1) Module 1" << endl;
    cout << "  (2) End program" << endl;
    
    do
    {
        cout << endl
             << "Please select the module you want to use: ";
        cin >> module;
    } while (module < 1 || module > 2);
    
    cout << endl; 
    cout << "Module selected = " << module << endl;

} // getModule function

// getOutMct
void getOutMct(int &outMct)
{
    do
    {
        cout << endl
             << "Write .mct output? (0 = no, 1 = yes): ";
        cin >> outMct;
    } while (outMct < 0 || outMct > 1);
	
} // getOutMct

// getOutMca
void getOutMca(int &outMca)
{
    do
    {
        cout << endl
             << "Write .mca output? (0 = no, 1 = yes): ";
        cin >> outMca;
    } while (outMca < 0 || outMca > 1);
	
} // getOutMca

// getOutMgd
void getOutMgd(int &outMgd)
{
    do
    {
        cout << endl
             << "Write .mgd output? (0 = no, 1 = yes): ";
        cin >> outMgd;
    } while (outMgd < 0 || outMgd > 1);
	
} // getOutMgd

// 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;
                    exit(0);                     
                }
                else
                {
                    if (countR <= dispProb[countD])
                    {
                        dispTable[countR] = dispStep[countD];
                        flagFoundD = 1;
                    }
                }
            
            } while(flagFoundD != 1);
        }
    }
    
    
} // getDispTable

// getNStep
int getNStep(int *dispTable)
{
    //
    int randInt;    // random int between 0  and maxRandomInt
    int nStepMoved; // number of steps moved
    
    randInt = getRandomInt(0, maxRandomInt);
    nStepMoved = dispTable[randInt];

    return nStepMoved;
} // getNStep


// checkParentPosition
void checkParentPosition(int &parentPosition, int optBoundary, int leftBoundary, int rightBoundary)
{
    //
	int flagWithinLimit = 0; 
	
	//
	do
	{
	    // if outside of the boundaries
		if (   parentPosition < leftBoundary || parentPosition > rightBoundary)
		{
	        switch(optBoundary)
	        {
	            // optBoundary == 1
	            // no boundaries, donut-shape habitat
	            case 1:     if (parentPosition < leftBoundary)
	                        {
	                            parentPosition = ((rightBoundary + 1) - (leftBoundary - parentPosition));
	                        }
	                        else
	                        {
	                            parentPosition = ((leftBoundary - 1) + (parentPosition - rightBoundary));
	                        }
	                        break;
	        
	            // optBoundary == 2
	            // absorbing boundaries
	            case 2:     if (parentPosition < leftBoundary)
	                        {
	                            parentPosition = leftBoundary;
	                        }
	                        else
	                        {
	                            parentPosition = rightBoundary;
	                        }
	                        break;

	            // optBoundary == 3
	            // reflecting boundaries
	            case 3:     if (parentPosition < leftBoundary)
	                        {
	                            parentPosition = ((leftBoundary - parentPosition) + leftBoundary);
	                        }
	                        else
	                        {
	                            parentPosition = (rightBoundary - (parentPosition - rightBoundary));
	                        }
	                        break;
	        } // end switch
		} // if outside of the boundaries
		else // within the boundaries
		{
			// reset flag
			flagWithinLimit = 1;
		}
	} while (flagWithinLimit == 0);
	
} // checkParentPosition


// getParentPosition
int getParentPosition(int countG, int selfPosition, int nPopSize, int optBoundary, int posBarrier, int durBarrier, int *dispTable)
{
	//
	int parentPosition;
    int nStepMoved;
	
    // get the number of steps moved
    nStepMoved = getNStep(dispTable);
    // calculate parent position
	parentPosition = selfPosition + nStepMoved;

	// check parent position
	if (countG > durBarrier) // no barrier
	{
	    // left boundary = 0
	    // right boundary = (nPopSize - 1)
	    checkParentPosition(parentPosition, optBoundary, 0, (nPopSize - 1));
	} // no barrier
	else // barrier exists
	{
	    if (selfPosition < posBarrier) // currrent position is in the left part
	    {
	        // left boundary = 0
	        // right boundart = (posBarrier - 1)
    	    checkParentPosition(parentPosition, optBoundary, 0, (posBarrier - 1));
	    }
	    else // current position is in the right part
	    {
	        // left boundary = posBarrier
	        // right boundart = (nPopSize - 1)
    	    checkParentPosition(parentPosition, optBoundary, posBarrier, (nPopSize - 1));
	    }
	} // barrier exists
		
	return parentPosition;
	
} // getParentPosition

// getCladeID
void getCladeID(int nSample, int countG, int *cladeID, int **coalTime)
{
	//
	int countS;
	
	// 
	for (countS = 0; countS < nSample; countS++)
	{
		if (coalTime[0][countS] < countG) // same clade as ind 0
		{
			cladeID[countS] = 0;
		}
		else // diff clade with ind 0
		{
			cladeID[countS] = 1;
		}
	}

} // getCladeID

// getClade0Count
void getClade0Count(int nLoc, int nInd, int *cladeID, int *clade0Count)
{
	//
	int countL;
	int countI;
	
	// reset
	for (countL = 0; countL < nLoc; countL++)
	{
		clade0Count[countL] = 0;	
	}
	
	
	// count
	for (countL = 0; countL < nLoc; countL++)
	{
		for (countI = 0; countI < nInd; countI++)
		{
			if (cladeID[(countL * nInd + countI)] == 0)
			{
				clade0Count[countL]++;
			}
			
		}	
	
	}
	

} // getClade0Count

// getUnsignCOA
float getUnsignCOA(int a, int b, int c, int d)
{
	//
	float unsignCOA;
	
	// 
	if (a == b) // if a==b, COA = 0
	{
	    unsignCOA = 0;
	}
	else // calculate unsignCOA
	{
    	unsignCOA = sqrt(float((a*d - b*c)*(a*d - b*c))/((a+b)*(c+d)*(a+c)*(b+d)));
	}

	return unsignCOA;

} // getUnsignCOA


// getArrayCOA
void getArrayCOA(int nLoc, int nInd, int *clade0Count, float **arrayCOA)
{
	//
	int countRow;
	int countCol;
	
	// 
	for (countRow = 0; countRow < (nLoc - 1); countRow++)
	{
		for (countCol = (countRow + 1); countCol < nLoc; countCol++)
		{
			arrayCOA[countRow][countCol] = getUnsignCOA(clade0Count[countRow], clade0Count[countCol], (nInd - clade0Count[countRow]), (nInd - clade0Count[countCol]));
		}
	
	}

	// set diagonal to 0
	for (countRow = 0; countRow < nLoc; countRow++)
	{
		arrayCOA[countRow][countRow] = 0;
	}

	//
	for (countRow = 0; countRow < (nLoc - 1); countRow++)
	{
		for (countCol = (countRow + 1); countCol < nLoc; countCol++)
		{
			arrayCOA[countCol][countRow] = arrayCOA[countRow][countCol];
		}
	
	}

	
} // getArrayCOA

// getAvCoalTime
float getAvCoalTime(int locI, int locJ, int nInd, int **coalTime)
{
	//
	int sumCoalTime = 0;
	float avCoalTime = 0;

	//
	int countRow;
	int countCol;
	
	// 
	for (countRow = 0; countRow < nInd; countRow++)
	{
		for (countCol = 0; countCol < nInd; countCol++)
		{
			sumCoalTime += coalTime[(locI * nInd + countRow)][(locJ * nInd + countCol)];
		}
	}

	avCoalTime = float (sumCoalTime) / (nInd * nInd);
	
	return avCoalTime;

} // getAvCoalTime

// getAvGD
void getAvGD(int nLoc, int nInd, int **coalTime, float **avGD)
{
	//
	int countRow;
	int countCol;
	
	// 
	for (countRow = 0; countRow < nLoc; countRow++)
	{
		for (countCol = countRow; countCol < nLoc; countCol++)
		{
			avGD[countRow][countCol] = getAvCoalTime(countRow, countCol, nInd, coalTime);
		}
	
	}

	//
	for (countRow = 0; countRow < (nLoc - 1); countRow++)
	{
		for (countCol = (countRow + 1); countCol < nLoc; countCol++)
		{
			avGD[countCol][countRow] = avGD[countRow][countCol];
		}
	
	}
} // getAvGD


// module1
void module1(void)
{

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

        // output file - log output 
            // declare fileName
            fileName outputNameLog; 
            // declare output file
            ofstream outputFileLog;
            // call openOutputLog function
            openOutputLog(outputFileLog, outputNameLog);
        // output file - Mct output 
            // declare fileName
            fileName outputNameMct; 
            // declare output file
            ofstream outputFileMct;
            // call openOutputMct function
            openOutputMct(outputFileMct, outputNameMct);
        // output file - Mca output 
            // declare fileName
            fileName outputNameMca; 
            // declare output file
            ofstream outputFileMca;
            // call openOutputMca function
            openOutputMca(outputFileMca, outputNameMca);
        // output file - Mgd output 
            // declare fileName
            fileName outputNameMgd; 
            // declare output file
            ofstream outputFileMgd;
            // call openOutputMgd function
            openOutputMgd(outputFileMgd, outputNameMgd);
        //
        writeOutputHeaderM1(outputFileLog, outputNameLog);

    // declare variables 
        // input settings
        int outMct = 0; // write Mct output
        int outMca = 0; // write Mca output
        int outMgd = 0; // write Mgd output
        
        // dsp input settings
        int nIter; // number of iterations to run
        int nPopSize; // population size
        int nLoc; // number of location sampled
        int nInd; // number of individuals sampled per location
        int optBoundary; // boundary option
        int posBarrier; // barrier position
        int durBarrier; // barrier duration
        int nDispStep;  // number of possible dispersal steps 
                                
        // variable
        int nRepInterval; // report interval
        int nSample; // total number of individuals sampled
        int distLoc; // distance between locations sampled
        int selfPosition; // hold the position of self
        int parentPosition; // hold the position of the parent
        
        //
        int flagRun;
        
        // loop counting
        int countI; // count Iteration
        int countG = 0; // count Generation
        int countL; // count location
        int countInd; // count individual
        int countS; // count Sample
        int countS1; // count Sample
        int countS2; // count Sample
        
    // get settings
    getOutMct(outMct);
    getOutMca(outMca);
    getOutMgd(outMgd);
   
    // readInputSettingDsp
    readInputSettingDsp(inputFileDsp, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);

    // calculate internal variables
    nRepInterval = nPopSize / 10; // report interval = N/10
    nSample=nLoc*nInd; // total number of individuals sampled
    distLoc=(nPopSize - nSample)/(nLoc - 1); // distance between locations sampled

    // write settings to log output file
    // call writeM1Setting function
    writeM1Setting(outputFileLog, outputNameLog, outputNameMct, outputNameMca, outputNameMgd, inputNameDsp, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);

    // declare and initiate array        
        // 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();
            
        // array to hold current position
        int *curPosition;
        
        curPosition = new int [nSample];
        if (curPosition == NULL)
            memoryError();
        
        // array to hold parent position
        int *parPosition;
        
        parPosition = new int [nSample];
        if (parPosition == NULL)
            memoryError();

        // array to hold id of individuals with shared parent
        int *idIndShared;
        
        idIndShared = new int [nSample];
        if (idIndShared == NULL)
            memoryError();

        // array to hold clade id of individuals 
        int *cladeID;
        
        cladeID = new int [nSample];
        if (cladeID == NULL)
            memoryError();

        // array to hold number of individuals in clade 0 for each location 
        int *clade0Count;
        
        clade0Count = new int [nLoc];
        if (clade0Count == NULL)
            memoryError();

        // array to hold coalescent time
        int **coalTime;
        
        coalTime = new int* [nSample];
        if (coalTime == NULL)
            memoryError();

        for (countS = 0; countS < nSample; countS++)
        {
            coalTime[countS] = new int [nSample];
            if (coalTime == NULL)
                memoryError();
        }        

        // array to hold coefficient of association
        float **arrayCOA;
        
        arrayCOA = new float* [nLoc];
        if (arrayCOA == NULL)
            memoryError();

		for (countL = 0; countL < nLoc; countL++)
		{
			arrayCOA[countL] = new float [nLoc];
	        if (arrayCOA == NULL)
	            memoryError();			 
		}
        
        // array to hold average genetic distance within/between locations
        float **avGD;
        
        avGD = new float* [nLoc];
        if (avGD == NULL)
            memoryError();

		for (countL = 0; countL < nLoc; countL++)
		{
			avGD[countL] = new float [nLoc];
	        if (avGD == NULL)
	            memoryError();			 
		}
        
            
       
    // read input file
    readInputDsp(inputFileDsp, nDispStep, dispStep, dispProb);

	// dispInputSettingDsp
    dispInputSettingDsp(inputNameDsp, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep, dispStep, dispProb);

	// writeOutputSetting
    writeOutputSetting(outputFileMct, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);
    writeOutputSetting(outputFileMca, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);
    writeOutputSetting(outputFileMgd, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);

    // getDispTable
    getDispTable(nDispStep, dispStep, dispProb, dispTable);
    
    // iteration loop
    for (countI = 0; countI < nIter; countI++)
    {
        // initiate data array values
            // curPosition
                countS=0;
                // 1 to (nLoc-1) locations
                for (countL = 0; countL < (nLoc - 1); countL++)
                {
                    for (countInd = 0; countInd < nInd; countInd++)
                    {
                        curPosition[countS] = countL*(nInd+distLoc) + countInd;
                        countS++;
                    }
                }
                // last location
                for (countInd = nInd; countInd > 0; countInd--)
                {
                    curPosition[countS] = (nPopSize - countInd);
                    countS++;
                }
                
            
            // coalTime
            for (countS1 = 0; countS1 < nSample; countS1++)
            {
                for (countS2 = countS1; countS2 < nSample; countS2++)
                {
                    if (countS1 == countS2) // coal time with self
                    {
                        coalTime[countS1][countS2]=0;
                    }
                    else // coal time with others
                    {
                        coalTime[countS1][countS2]=1;
                    }
                }
            
            }
        
        // reset countG
        countG=0;
        
        // generation loop, go backward in time
        do
        {
            // reset flag
            flagRun = 0;
            
            
            // reset idIndShared
            for (countS = 0; countS < nSample; countS++)
            {
            	idIndShared[countS] = -1;
            }
            
            // check for individuals with shared parent
            for (countS1 = (nSample - 1); countS1 > 0; countS1--)
            {
	            for (countS2 = 0; countS2 < countS1; countS2++)
	            {
	            	if (curPosition[countS1] == curPosition[countS2])
	            	{
	            		idIndShared[countS1] = countS2;
	            		break;
	            	}
	            
	            }
            
            }
            
            // find parPosition
            for (countS = 0; countS < nSample; countS++)
            {
            	// check if shared parent with others
            	if (idIndShared[countS] == -1) // did not share parent with others
            	{
            		// find parPosition
            		selfPosition = curPosition[countS];
                    parentPosition = getParentPosition(countG, selfPosition, nPopSize, optBoundary, posBarrier, durBarrier, dispTable);
					parPosition[countS] = parentPosition;
            	}
            	else // shared parent with others
            	{
            		parPosition[countS] = parPosition[(idIndShared[countS])];
            	}
            } // find parPosition
            
            // update coalTime
            for (countS1 = 0; countS1 < (nSample - 1); countS1++)
            {
                for (countS2 = (countS1 + 1); countS2 < nSample; countS2++)
                {
                    if (parPosition[countS1] != parPosition[countS2]) // do not share parent
                    {
                        coalTime[countS1][countS2]++;
                        flagRun = 1;
                    }
                }            
            }
   
			// update curPosition for next loop
            for (countS = 0; countS < nSample; countS++)
            {
            	curPosition[countS] = parPosition[countS];
            }
			            
            
            // update countG (generation count)
            countG++;
            
            // report to screen
            if ((countG % nRepInterval) == 0)
            {
            	cout << "countI = " << countI 
            		 << ", countG = " << countG
            		 << endl;
            }
        } while (flagRun != 0);// end of generation loop
    
    	// fill coalTime array (below diag)
        for (countS1 = 0; countS1 < (nSample-1); countS1++)
        {
            for (countS2 = (countS1+1); countS2 < nSample; countS2++)
            {
            	coalTime[countS2][countS1] = coalTime[countS1][countS2];
            }        
        }
        
        // getCladeID
        getCladeID(nSample, countG, cladeID, coalTime);
        
        // getClade0Count
        getClade0Count(nLoc, nInd, cladeID, clade0Count);
        
		// getArrayCOA
		getArrayCOA(nLoc, nInd, clade0Count, arrayCOA);

		// getAvGD
		getAvGD(nLoc, nInd, coalTime, avGD);

        // write output
        // Mct
        if (outMct == 1)
        {
            writeMct(outputFileMct, countI, countG, nSample, coalTime);
        }

        // Mca
        if (outMca == 1)
        {
            writeMca(outputFileMca, countI, countG, nLoc, arrayCOA);
        }
        
        // Mgd
        if (outMgd == 1)
        {
            writeMgd(outputFileMgd, countI, countG, nLoc, avGD);
        }
        
    
    } // iteration loop


	// delete data arrays
	// release memory
	delete [] dispStep;
	delete [] dispProb;
	delete [] dispTable;
	delete [] curPosition;
	delete [] parPosition;
	delete [] idIndShared;
	delete [] cladeID;
	delete [] clade0Count;
	delete [] coalTime;
	delete [] arrayCOA;
	delete [] avGD;
        
    
    // close input/output files
    // call closeFile function
        // close input file
        closeInputFile(inputFileDsp);
        // close output file
        closeOutputFile(outputFileLog);
        closeOutputFile(outputFileMct);
        closeOutputFile(outputFileMca);

} // module1