// breakAna 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] = "breakAna";
    // program name
const char versionNumber[15] = "A2";
    // version number
const char versionDate[15] = "04/15/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 = 10000;
    // 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


// openInputFileMgd function
// open input file and get input file name
void openInputFileMgd(ifstream &inputFileMgd, fileName &inputNameMgd)
{
 
    // 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 (Genetic distance matrix <.mgd>): ";
    cin >> inputNameMgd;

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

} // end of openInputFileMgd


// 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 (breakAna Log <.bal>): ";
    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


// 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();
}


// readInputSetting function
void readInputSetting(ifstream &inputFile, 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
    {
        inputFile.get(symbolInput);
        switch (symbolInput)
        {
            case symbolComment:     inputFile.getline(inputComment, 81, '\n');
                                    break;
            case symbolNIter:       inputFile >> nIter;
                                    flagNIter = 1;
                                    break;
            case symbolOptSize:     inputFile >> nPopSize;
                                    inputFile >> nLoc;
                                    inputFile >> nInd;                                    
                                    flagOptSize = 1;
                                    break;
            case symbolOptBoundary: inputFile >> optBoundary;
                                    inputFile >> posBarrier;
                                    inputFile >> durBarrier;
                                    flagOptBoundary = 1;
                                    break;
            case symbolNDispStep:   inputFile >> 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);
    }
    
} // readInputSetting function


// readInputMgd
void readInputMgd(ifstream &inputFile, int countI, int nLoc, int *maxCoalTime, float **avGD)
{    
    // declare variables
    char symbolInput;
    int iterID;
    char inputComment[81];
    int flagGotGD = 0;
    
    //
    int countRow;
    int countCol;
    
    // read input setting
    do
    {
        inputFile.get(symbolInput);
        switch (symbolInput)
        {
            case symbolComment:     inputFile.getline(inputComment, 81, '\n');
                                    break;
            case symbolMGD:         inputFile >> iterID;
                                    inputFile >> maxCoalTime[countI];
                                    for (countRow = 0; countRow < nLoc; countRow++)
                                    {
                                        for (countCol = 0; countCol < nLoc; countCol++)
                                        {
                                            inputFile >> avGD[countRow][countCol];
                                        }                                    
                                    }
                                    flagGotGD = 1;
                                    break;
            default:                flagGotGD = 0;
                                    break;
        }
    } while (flagGotGD !=1 && !inputFile.eof());
    //

} // readInputMgd

// dispInputSetting
void dispInputSetting(fileName &inputName, int nIter, int nPopSize, int nLoc, int nInd, int optBoundary, int posBarrier, int durBarrier, int nDispStep)
{
	
	//
	cout << endl;
	cout << "Input setting from " << inputName << 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;
	
} // dispInputSetting


// 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 inputNameMgd, 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 << "   Mgd file: " << inputNameMgd << endl;
    outputFileLog << " Output file" << endl;
    outputFileLog << "   Log file: " << outputNameLog << 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


// writeBreakLoc
void writeBreakLoc(ofstream &outputFile, fileName inputNameMgd, int nLoc, int *breakLocCount, float *breakLocFreq)
{
	//
	int countL;
	
	//
	outputFile << "*** Start of BreakLoc (" << inputNameMgd << ")***" << endl;
	
	outputFile << inputNameMgd << "\t" << "Loc" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << countL << "\t";
	}
	outputFile << endl;
	
	outputFile << inputNameMgd << "\t" << "Count" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << breakLocCount[countL] << "\t";
	}
	outputFile << endl;
	
	outputFile << inputNameMgd << "\t" << "Freq" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << setprecision(4)
                   << setiosflags(ios::fixed | ios::showpoint)
                   << breakLocFreq[countL]
                   << "\t";
	}
	outputFile << endl;
	
	outputFile << "*** End of BreakLoc (" << inputNameMgd << ")***" << endl;
	outputFile << endl;

} // writeBreakLoc

// writeBreakDist
void writeBreakDist(ofstream &outputFile, fileName inputNameMgd, int nLoc, int *breakDistCount, float *breakDistFreq)
{
	//
	int countL;
	
	//
	outputFile << "*** Start of BreakDist (" << inputNameMgd << ")***" << endl;
	
	outputFile << inputNameMgd << "\t" << "Dist" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << countL << "\t";
	}
	outputFile << endl;
	
	outputFile << inputNameMgd << "\t" << "Count" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << breakDistCount[countL] << "\t";
	}
	outputFile << endl;
	
	outputFile << inputNameMgd << "\t" << "Freq" << "\t";
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
	    outputFile << setprecision(4)
                   << setiosflags(ios::fixed | ios::showpoint)
                   << breakDistFreq[countL]
                   << "\t";
	}
	outputFile << endl;
	
	outputFile << "*** End of BreakDist (" << inputNameMgd << ")***" << endl;
	outputFile << endl;

} // writeBreakDist


// 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



// getAvGDNeiLoc
void getAvGDNeiLoc(int countI, int nLoc, float **avGD, float **avGDNeiLoc)
{
	// 
	int countL;
	
	//
	for (countL = 0; countL < (nLoc - 1); countL++)
	{
		avGDNeiLoc[countI][countL] = avGD[countL][(countL + 1)];
	}

} // getAvGDNeiLoc


// getMaxBreak
void getMaxBreak(int nIter, int nLoc, float **avGDNeiLoc, float *maxBreakStr, int *maxBreakLoc)
{
	//
	float maxStr;
	int maxLoc;
	
	//
	int countI;
	int countL;
	
	//
	for (countI = 0; countI < nIter; countI++)
	{
		// reset
		maxStr = 0;
		maxLoc = 0;
		
		// find maxBreak
		for (countL = 0; countL < (nLoc - 1); countL++)
		{
			if (avGDNeiLoc[countI][countL] > maxStr)
			{
				maxStr = avGDNeiLoc[countI][countL];
				maxLoc = countL;
			}
		}
		
		maxBreakStr[countI] = maxStr;
		maxBreakLoc[countI] = maxLoc;
	}
	

} // getMaxBreak


// getBreakLoc
void getBreakLoc(int nIter, int nLoc, int *maxBreakLoc, int *breakLocCount, float *breakLocFreq)
{
    //
    int countI;
    int countL;
    
    //
    for (countI = 0; countI < nIter; countI++)
    {
        breakLocCount[(maxBreakLoc[countI])]++;
    }
    
    //
    for (countL = 0; countL < (nLoc - 1); countL++)
    {
        breakLocFreq[countL] = float(breakLocCount[countL]) / nIter;
    }


} // getBreakLoc

// getBreakDist
void getBreakDist(int nIter, int nLoc, int *maxBreakLoc, int *breakDistCount, float *breakDistFreq)
{
    //
    int nPair = nIter * (nIter + 1) / 2; // number of tree pairs 
    int breakLoc1;
    int breakLoc2;
    int breakLocDist;
    
    //
    int countI1;
    int countI2;
    int countL;
    
    //
    for (countI1 = 0; countI1 < nIter; countI1++)
    {
        for (countI2 = countI1; countI2 < nIter; countI2++)
        {
            //
            breakLoc1 = maxBreakLoc[countI1];
            breakLoc2 = maxBreakLoc[countI2];
    		
    		//
    		breakLocDist = abs(breakLoc1 - breakLoc2);
    		
    		//
			breakDistCount[breakLocDist]++;
        }    
    }

    //
    for (countL = 0; countL < (nLoc - 1); countL++)
    {
        breakDistFreq[countL] = float(breakDistCount[countL]) / nPair;
    }


} // getBreakDist



// module1
void module1(void)
{

    // file operation
        // input file - Mgd
            // declare fileName
            fileName inputNameMgd; 
            // declare input file
            ifstream inputFileMgd;
            // call openInputFile function
            openInputFileMgd(inputFileMgd, inputNameMgd);

        // output file - log output 
            // declare fileName
            fileName outputNameLog; 
            // declare output file
            ofstream outputFileLog;
            // call openOutputLog function
            openOutputLog(outputFileLog, outputNameLog);
        //
        writeOutputHeaderM1(outputFileLog, outputNameLog);

    // declare variables 
        // user input settings
        
        // input file 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
        //
        
        // loop counting
        int countI; // count Iteration
        int countL; // count location
        
   
    // readInputSetting
    readInputSetting(inputFileMgd, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);


	// dispInputSetting
    dispInputSetting(inputNameMgd, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);
    
    // calculate internal variables

    // write settings to log output file
    // call writeM1Setting function
    writeM1Setting(outputFileLog, outputNameLog, inputNameMgd, nIter, nPopSize, nLoc, nInd, optBoundary, posBarrier, durBarrier, nDispStep);
    
    // declare and initiate array        

        // array to hold countG (max coal time)
        int *maxCoalTime;

        maxCoalTime = new int [nIter];
        if (maxCoalTime == NULL)
            memoryError();
            
        // array to hold avGD array from input file
        // [nLoc][nLoc]
        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();			 
		}

        // array to hold avGD between neighbor locations
        // [nIter][(nLoc - 1)]
        float **avGDNeiLoc;
        
        avGDNeiLoc = new float* [nIter];
        if (avGDNeiLoc == NULL)
            memoryError();

		for (countI = 0; countI < nIter; countI++)
		{
			avGDNeiLoc[countI] = new float [(nLoc - 1)];
	        if (avGDNeiLoc == NULL)
	            memoryError();			 
		}
		
		// array to hold max break strength of all input
		float *maxBreakStr;
		
        maxBreakStr = new float [nIter];
        if (maxBreakStr == NULL)
            memoryError();
		
		// array to hold max break location of all input
		int *maxBreakLoc;
		
        maxBreakLoc = new int [nIter];
        if (maxBreakLoc == NULL)
            memoryError();
		
		// array to hold count of break location
		int *breakLocCount;
		
        breakLocCount = new int [(nLoc - 1)];
        if (breakLocCount == NULL)
            memoryError();
            
            // reset
			for (countL = 0; countL < (nLoc - 1); countL++)
			{
				breakLocCount[countL] = 0;
			}
            
		// array to hold freq of break location
		float *breakLocFreq;
		
        breakLocFreq = new float [(nLoc - 1)];
        if (breakLocFreq == NULL)
            memoryError();
                                
		// array to hold count of break distance
		int *breakDistCount;
		
        breakDistCount = new int [(nLoc - 1)];
        if (breakDistCount == NULL)
            memoryError();
            
            // reset
			for (countL = 0; countL < (nLoc - 1); countL++)
			{
				breakDistCount[countL] = 0;
			}
            
		// array to hold freq of break distance
		float *breakDistFreq;
		
        breakDistFreq = new float [(nLoc - 1)];
        if (breakDistFreq == NULL)
            memoryError();
                                
       
    // read input file
    for (countI = 0; countI < nIter; countI++)
    {
		readInputMgd(inputFileMgd, countI, nLoc, maxCoalTime, avGD);

		getAvGDNeiLoc(countI, nLoc, avGD, avGDNeiLoc);
    }
    
    // getMaxBreak
	getMaxBreak(nIter, nLoc, avGDNeiLoc, maxBreakStr, maxBreakLoc);
    
    // getBreakLoc
    getBreakLoc(nIter, nLoc, maxBreakLoc, breakLocCount, breakLocFreq);
    
    // getBreakDist
    getBreakDist(nIter, nLoc, maxBreakLoc, breakDistCount, breakDistFreq);

    // writeBreakLoc
    writeBreakLoc(outputFileLog, inputNameMgd, nLoc, breakLocCount, breakLocFreq);

    // writeBreakDist
    writeBreakDist(outputFileLog, inputNameMgd, nLoc, breakDistCount, breakDistFreq);
    
	// delete data arrays
	// release memory
	delete [] maxCoalTime;
	delete [] avGD;
	delete [] avGDNeiLoc;
	delete [] maxBreakStr;
	delete [] maxBreakLoc;
	delete [] breakLocCount;
	delete [] breakLocFreq;
	delete [] breakDistCount;
	delete [] breakDistFreq;
        
    
    // close input/output files
    // call closeFile function
        // close input file
        closeInputFile(inputFileMgd);
        // close output file
        closeOutputFile(outputFileLog);

} // module1