//  Modified version of Paul Lewis' ncltest program that will rearrange taxa in specific orders
//

#include <string.h>
#include "ncl.h"
#include "ncltest.h"

#define MAXSITE 200000

void usage() {
	cerr << "Usage: ./id_gappy <infile> <outfile> <gap/ambiguous_taxa> <length_flag>\n\n";
	cerr << "  For <gap/ambiguous_taxa> use:\n";
	cerr << "  - use <g=#> for the maximum number of taxa with gaps or ambiguous\n";
	cerr << "    states (i.e., columns will be placed in the gappy set if they have\n";
	cerr << "    more than two gap/ambiguous state elements if g=2).\n";
	cerr << "  - use <d=#> for the minimum number taxa with defined (non-gap, non-\n";
	cerr << "    ambiguous states (i.e., columns will be placed in the gappy set if\n";
	cerr << "    they have do not have at least two defined elements if d=2).\n\n";
	cerr << "  <length_flag> refers to the length of 'gappy' regions considered long.\n";
	cerr << "  Long gappy regions will be listed in the terminal window. Note that\n";
	cerr << "  the <length_flag> is optional, and it should be an integer.\n\n";
}

/*----------------------------------------------------------------------------------------------------------------------
|	Returns a true if the character passed is a number
*/
bool isnumeric(char iv) {
	bool num=false;
	if (iv=='0' || iv=='1') num=true;
	if (iv=='2' || iv=='3') num=true;
	if (iv=='4' || iv=='5') num=true;
	if (iv=='6' || iv=='7') num=true;
	if (iv=='8' || iv=='9') num=true;
	return num;
}

/*----------------------------------------------------------------------------------------------------------------------
|	Returns a true if the character passed is a number
*/
int char2int(char iv) {
	int v;
	if (iv=='0') v=0;
	if (iv=='1') v=1;
	if (iv=='2') v=2;
	if (iv=='3') v=3;
	if (iv=='4') v=4;
	if (iv=='5') v=5;
	if (iv=='6') v=6;
	if (iv=='7') v=7;
	if (iv=='8') v=8;
	if (iv=='9') v=9;
	return v;
}

/*----------------------------------------------------------------------------------------------------------------------
|	Returns a positive integer if the string passed beings with "g=" or "d=", otherwise returns -1
*/
int getgapparam(char* in) {
	int val=-1;
	bool problem=false;
	unsigned inlen=strlen(in);
	if (in[0]=='g' || in[0]=='d') {
		if (inlen>2 && isnumeric(in[2])) {
			val=char2int(in[2]);
		} // end if (inlen>2...
		else problem=true;
		if (inlen>3) {
			for (unsigned ggp=3; ggp<inlen; ggp++) {
				if (isnumeric(in[ggp])) {
					val=val*10;
					val=val+char2int(in[ggp]);
				}
				else {
					problem=true;
					ggp=inlen;
				}
			} // end for (unsigned ggp=2; ggp<inlen...
		} // end if (inlen>3...
	} // end if (in[0]=='g' || in[0]=='d'...
	if (problem) val=-1;
	return val;
}

/*----------------------------------------------------------------------------------------------------------------------
|	Returns an IUPAC ambiguity code representing the char* in (which will contain A, C, G, and T) for sites
|   based upon the possible nucleotides found at a specific site
*/
char ambigDNA(char* in, unsigned s) {
  char st[256];
  strcpy(st,in);
  char val;
  if (s==2) {
    if (st[0]=='A') {
	  if (st[1]=='C') val='M';
	  else if (st[1]=='G') val='R';
		   else val='W'; // must be A and T
	}
	else if (st[0]=='C') {
	       if (st[1]=='G') val='S';
		   else val='Y'; // must be C and T
	     }
		 else val='K';  // must be G and T
  }
  else if (s==3) {
         if (st[0]=='A') {
		   if (st[1]=='C') {
		     if (st[3]=='G') val='V'; // data is A, C, or G
			 else val='H'; // data must be A, C, or T
		   }
		   else val='D'; // data must be A, G, or T
		 }
		 else val='B'; // data must be C, G, or T
       }
       else {
	     val='N'; // data must be A, C, G, or T
	   }
  return val;
}

/*----------------------------------------------------------------------------------------------------------------------
|	Takes an ostream reference in addition to the istream reference required by the base class NxsToken. The istream 
|	'is' is passed along to the base class constructor, and 'os' is used to initialize the data member out.
*/
MyToken::MyToken(
  istream &is,	/* reference to input file stream attached to the NEXUS data file to be read */
  ostream &os)	/* reference to the output stream used for outputting results */
  : out(os), NxsToken(is)
	{
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Overrides the base class version. When a NEXUS printable comment is encountered in a NEXUS file, the comment is 
|	passed to the virtual function OutputComment in its argument 'msg'. This version outputs 'msg' directly to both the
|	standard output (cout) and the output file (out).
*/
void MyToken::OutputComment(
  const NxsString &msg)	/* reference to string containing the printable comment embedded in the NEXUS data file */
	{
	cout << msg << endl;
	out << msg << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Initializes infile_exists to true if a file by the name 'infname' exists and to false otherwise. If a file named
|	'infname' exists, proceeds to open the ifstream data member inf to that file and the ofstream data member outf to
|	the file specified by the 'outfname' argument. Does not check for the existance of the output file beforehand, so
|	will overwrite a file if one by than name already exists. In a real application, use the FileExists member function
|	to check for the existance of the output file as well and warn the user before opening it.
*/
MyNexusFileReader::MyNexusFileReader(
  char *infname,	/* the name of the input file (should be a NEXUS-formatted data file) */
  char *outfname)	/* the name of a file to send output to */
 : NxsReader() 
	{
	infile_exists = true;
	if (!FileExists(infname))
		infile_exists = false;

	if (infile_exists)
		{
		inf.open(infname, ios::in | ios::binary);
		outf.open(outfname);
		}
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Destructor simply closes the input file stream inf and the output file stream outf.
*/
MyNexusFileReader::~MyNexusFileReader() 
	{
	inf.close();
	outf.close();
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports to both standard output and the output file stream outf that
|	a block named 'blockName' was just entered by the NEXUS file reader. Always returns true.
*/
bool MyNexusFileReader::EnteringBlock(
  NxsString blockName)	/* the name of the block just entered */
	{
	cout << "Reading \"" << blockName << "\" block..." << endl;
	outf << "Reading \"" << blockName << "\" block..." << endl;
	return true;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports to both standard output and the output file stream outf that
|	the NEXUS file reader just finished with a block named 'blockName'.
*/
void MyNexusFileReader::ExitingBlock(
  NxsString blockName)	/* the name of the block being exited */
	{
	cout << "Finished with \"" << blockName << "\" block." << endl;
	outf << "Finished with \"" << blockName << "\" block." << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply outputs the string 's' to both the standard output and the output
|	stream outf. This function will be called when a printable NEXUS comment is found in a NEXUS file. The contents of
|	the comment are provided in the string argument 's'.
*/
void MyNexusFileReader::OutputComment(
  const NxsString &s)	/* the contents of the printable NEXUS comment */
	{
	cout << endl;
	cout << s << endl;
	outf << endl;
	outf << s << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	A simple function to test for the existance of a file named 'fn'. Returns true if an attempt to open the file 
|	succeeds, and false otherwise. No attempt is made to diagnose failure to open the file, so this function may return
|	false even when the file exists if the file could not be opened for some other reason (e.g. user does not have
|	read permissions for that file).
*/
bool MyNexusFileReader::FileExists(
  char *fn)	/* the name of the file to test (path should be included if you do not want to assume the current working directory) */
	{
	if (fn == NULL) 
		return false;

	bool exists = false;

	FILE *f = fopen(fn, "r");
	if (f != NULL)
		{
		exists = true;
		fclose(f);
		f = NULL;
		}

	return exists;
	}
 
/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports to both the standard output and to the output stream outf
|	that a block named 'blockName' is being skipped because this NEXUS file reader does not know how to interpret blocks
|	by this name. This situation is common if the NEXUS file has private blocks, e.g. a PAUP block containing commands
|	understood only by the program PAUP.
*/
void MyNexusFileReader::SkippingBlock(
  NxsString blockName)	/* the name of the block being skipped */
	{
	cout << "Skipping unknown block (" << blockName << ")..." << endl;
	outf << "Skipping unknown block (" << blockName << ")..." << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports to both the standard output and to the output stream outf
|	that a block named 'blockName' is being skipped because this block, though known to the NEXUS reader, is currently
|	disabled. This situation arises when reading a treefile for instance (the application may disable all blocks except
|	NxsTreesBlock because only TREES blocks expected in a treefile. Overriding this function gives the application a
|	chance to report to the user that NEXUS blocks were found in the file that were unexpected under the circumstances.
*/
void MyNexusFileReader::SkippingDisabledBlock(
  NxsString blockName)	/* the name of the disabled block being skipped */
	{
	cout << "Skipping disabled block (" << blockName << ")..." << endl;
	outf << "Skipping disabled block (" << blockName << ")..." << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. If a block of the type referenced by 'nexusBlock' has been read (i.e. the
|	IsEmpty member function returns false), this function calls the Report member function of 'nexusBlock' to dump the
|	contents of the block object to the output file stream outf.
*/
void MyNexusFileReader::DebugReportBlock(
  NxsBlock &nexusBlock)	/* reference to the NxsBlock to be dumped if not empty */
	{
	if (!nexusBlock.IsEmpty())
		{
		outf << endl;
		outf << "********** Contents of the " << nexusBlock.GetID() << " block **********" << endl;
		nexusBlock.Report(outf);
		}
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. This function simply reports the error message contained in the string
|	'msg' and the file position, line and column (contained in 'pos', 'line' and 'col', respectively) to both the 
|	standard error stream and the output file stream outf. It then aborts the application by invoking exit(0).
*/
void MyNexusFileReader::NexusError(
  NxsString msg,	/* reference to the error message string */
  file_pos pos,		/* the file position where the error was discovered */
  unsigned line,		/* the line where the error was discovered */
  unsigned col)			/* the column where the error was discovered */
	{
	cerr << endl;
	cerr << "Error found at line " << line;
	cerr << ", column " << col;
	cerr << " (file position " << pos << "):" << endl;
	cerr << msg << endl;

	outf << endl;
	outf << "Error found at line " << line;
	outf << ", column " << col;
	outf << " (file position " << pos << "):" << endl;
	outf << msg << endl;

	cerr << endl;
	cerr << "Press return to quit..." << endl;
	cin.get();

	exit(0);
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Class derived from NxsCharactersBlock so that an output file stream can be used for reporting unknown commands that
|	are being skipped. Overrides SkippingCommand function for this purpose.
*/
MyCharactersBlock::MyCharactersBlock(
  NxsTaxaBlock *tb,			/* pointer to the NxsTaxaBlock object to use to store taxon labels */
  NxsAssumptionsBlock *ab,	/* pointer to the NxsAssumptionsBlock object to use for storing charsets, taxsets, etc. */
  ostream &o)				/* reference to the ostream object for use in reporting commands that are being skipped */
  : outf(o), NxsCharactersBlock(tb, ab)
	{
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports the command being skipped (name of command is in 's') both 
|	to the standard output stream as well as to the output stream outf.
*/
void MyCharactersBlock::SkippingCommand(
  NxsString s)	/* the name of the command being skipped */
	{
	cout << "Skipping unknown command (" << s << ")..." << endl;
	outf << "Skipping unknown command (" << s << ")..." << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Class derived from NxsDataBlock so that an output file stream can be used for reporting unknown commands that are 
|	being skipped. Overrides SkippingCommand function for this purpose.
*/
MyDataBlock::MyDataBlock(
  NxsTaxaBlock *tb,			/* pointer to the NxsTaxaBlock object to use to store taxon labels */
  NxsAssumptionsBlock *ab,	/* pointer to the NxsAssumptionsBlock object to use for storing charsets, taxsets, etc. */
  ostream &o)				/* reference to the ostream object for use in reporting commands that are being skipped */
  : outf(o), NxsDataBlock(tb, ab)
	{
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports the command being skipped (name of command is in 's') both 
|	to the standard output stream as well as to the output stream outf.
*/
void MyDataBlock::SkippingCommand(
  NxsString s)	/* the name of the command being skipped */
	{
	cout << "Skipping unknown command (" << s << ")..." << endl;
	outf << "Skipping unknown command (" << s << ")..." << endl;
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Class derived from NxsAssumptionsBlock so that an output file stream can be used for reporting unknown commands 
|	that are being skipped. Overrides SkippingCommand function for this purpose.
*/
MyAssumptionsBlock::MyAssumptionsBlock(
  NxsTaxaBlock *tb,			/* pointer to the NxsTaxaBlock object to use to store taxon labels */
  ostream &o)				/* reference to the ostream object for use in reporting commands that are being skipped */
  : outf(o), NxsAssumptionsBlock(tb)
	{
	}

/*----------------------------------------------------------------------------------------------------------------------
|	Override of base class virtual function. Simply reports the command being skipped (name of command is in 's') both 
|	to the standard output stream as well as to the output stream outf.
*/
void MyAssumptionsBlock::SkippingCommand(
  NxsString s)	/* the name of the command being skipped */
	{
	cout << "Skipping unknown command (" << s << ")..." << endl;
	outf << "Skipping unknown command (" << s << ")..." << endl;
	}

int main(int argc, char *argv[]) {

	if (argc<4 || argc>5) {
		usage();
		return 1;
	}
	int critlen=-1;
	if (argc==5) {
		critlen=atoi(argv[4]);
	}

	char infname[256];
	char outfname[256];
	char newfname[256];
	char authfname[256];
	strcpy(infname,argv[1]);
	strcpy(outfname,"nexus.log");
	strcpy(newfname,argv[2]);
	
	char gapparam[256];
	strcpy(gapparam,argv[3]);
	int maxgap=getgapparam(gapparam);
	bool gapdef=false;
	if (gapparam[0]=='g' && maxgap!=-1) {
		cerr << "Maximum numbers of gaps is ";
		gapdef=true;
	}
	if (gapparam[0]=='d' && maxgap!=-1) {
		cerr << "Minimum number of defined sites is ";
		gapdef=true;
	}
	if (!gapdef) {
		cerr << "Sites included in 'gappy' set are defined incorrectly.\n\n";
		usage();
		return 1;
	}
	else cerr << maxgap << ".\n\n";

	MyNexusFileReader nexus(infname, outfname);
	if (!nexus.infile_exists) {
		ofstream outf(outfname);
		outf << "Error: specified nexus input file (";
		outf << infname;
		outf << ") does not exist." << endl;
		outf.close();

		cerr << "Error: specified nexus input file (";
		cerr << infname;
		cerr << ") does not exist." << endl;

		cerr << endl;
		cerr << "Press return to quit..." << endl;
		cin.get();

		return 0;
	}

	NxsTaxaBlock		*taxa			= new NxsTaxaBlock();
	MyAssumptionsBlock	*assumptions	= new MyAssumptionsBlock(taxa, nexus.outf);
	NxsTreesBlock		*trees			= new NxsTreesBlock(taxa);
	MyCharactersBlock	*characters		= new MyCharactersBlock(taxa, assumptions, nexus.outf);
	MyDataBlock			*data			= new MyDataBlock(taxa, assumptions, nexus.outf);
	NxsDistancesBlock	*distances		= new NxsDistancesBlock(taxa);

	nexus.Add(taxa);
	nexus.Add(assumptions);
	nexus.Add(trees);
	nexus.Add(characters);
	nexus.Add(data);
	nexus.Add(distances);

	MyToken token(nexus.inf, nexus.outf);
	nexus.Execute(token);

	nexus.outf << endl;
	nexus.outf << "Reports follow on all blocks currently in memory." << endl;

	if (!data->IsEmpty()) {
		data->Report(nexus.outf);
		data->TransferTo(*characters);
	}
	if (!taxa->IsEmpty())
		taxa->Report(nexus.outf);
	if (!trees->IsEmpty())
		trees->Report(nexus.outf);
	if (!characters->IsEmpty())
		characters->Report(nexus.outf);
	if (!distances->IsEmpty())
		distances->Report(nexus.outf);
	if (!assumptions->IsEmpty())
		assumptions->Report(nexus.outf);

	NxsStringVector z;
	nexus.outf << endl;
	if (assumptions->GetNumTaxSets() > 0)
		{
		assumptions->GetTaxSetNames(z);
		NxsStringVector::const_iterator zi;
		for (zi = z.begin(); zi != z.end(); zi++)
			{
			nexus.outf << "Taxset " << (*zi) << ": ";
			NxsUnsignedSet &x = assumptions->GetTaxSet(*zi);
			NxsUnsignedSet::const_iterator xi;
			for (xi = x.begin(); xi != x.end(); xi++)
				nexus.outf << (*xi + 1) << " ";
			nexus.outf << endl;
			}
		}

	if (assumptions->GetNumCharSets() > 0)
		{
		assumptions->GetCharSetNames(z);
		NxsStringVector::const_iterator zi;
		for (zi = z.begin(); zi != z.end(); zi++)
			{
			nexus.outf << "Charset " << (*zi) << ": ";
			NxsUnsignedSet &x = assumptions->GetCharSet(*zi);
			NxsUnsignedSet::const_iterator xi;
			for (xi = x.begin(); xi != x.end(); xi++)
				nexus.outf << (*xi + 1) << " ";
			nexus.outf << endl;
			}
		}

	if (assumptions->GetNumExSets() > 0)
		{
		assumptions->GetExSetNames(z);
		NxsStringVector::const_iterator zi;
		for (zi = z.begin(); zi != z.end(); zi++)
			{
			nexus.outf << "Exset " << (*zi) << ": ";
			NxsUnsignedSet& x = assumptions->GetExSet(*zi);
			NxsUnsignedSet::const_iterator xi;
			for (xi = x.begin(); xi != x.end(); xi++)
				nexus.outf << (*xi + 1) << " ";
			nexus.outf << endl;
			}
		}

	ofstream logf("id_gappy.log");
	unsigned i, j, k, currChar, numgaps, gapstart, gappylen, output, mindef;
	int ntax = characters->GetNTax();
	cout << "Ntax  = " << ntax << endl;
	logf << "Ntax  = " << ntax << endl;
	int nchar = characters->GetNChar();
	cout << "Nchar = " << nchar << endl;
	logf << "Nchar = " << nchar << endl;
	if (gapparam[0]=='d') {
		mindef=maxgap;
		maxgap=ntax-mindef-1;
	}
	
	bool gapsites[MAXSITE];
	for (i=0; i<nchar; i++) gapsites[i]=false;
	/*----------------------------------------------------------------------------------------------------------------------
	|  Find the gappy sites and output a CHARSET.
	*/
	output=0;
	ofstream newf(newfname);
	newf << "\nBegin Sets;\n";
	newf << "[ Two versions of the same charset are presented below - delete one ]\n\n";
	if (gapparam[0]=='g') newf << "\tCHARSET gapsites_g" << maxgap << " = ";
	else newf << "\tCHARSET gapsites_d" << mindef << " = ";
	
	for (currChar = 0; currChar < nchar; currChar++) {
		numgaps=0;
		logf << currChar+1 << " - ";
		for (i = 0; i < ntax; i++) {
			if (characters->IsGapState(i, currChar)) {
				numgaps++;
				logf << '-';
			}
			else if (characters->IsMissingState(i, currChar)) {
					numgaps++;
					logf << '?';
				}
				else if (characters->GetNumStates(i, currChar)) logf << characters->GetState(i, currChar, 0);
					else logf << 'N';
		} // for (i = 0;...
		logf << " - " << numgaps;
	// 	cerr << currChar+1 << '-' << numgaps << ' ';  // - Direct output for debugging
		if (numgaps>maxgap) { 
			logf << "*";
			gapsites[currChar]=true;
			newf << currChar+1 << ' ';
			output++;
			if (output==10) {
				newf << endl << "\t\t\t";
				output=0;
			}
		}
		logf << endl;
	} // for (unsigned currChar = 0...
	newf << ";\n\n";
	logf.close();
	
	output=0;
	if (gapparam[0]=='g') newf << "\tCHARSET gapsites_g" << maxgap << " = ";
	else newf << "\tCHARSET gapsites_d" << mindef << " = ";
	for (unsigned gs=0; gs<nchar; gs++) {
		if (gapsites[gs]) {
			newf << gs+1 << ' ';
			if (gapsites[gs+1]) {
				gapstart=gs+1;
				gappylen=1;
				newf << "- ";
				while (gapsites[gs]) {
					gs++;
					gappylen++;
				}
				newf << gs << "  ";
				if (gappylen>critlen) cout << "  Long gap (" << gappylen << " sites) starts at POSITION " << gapstart << endl;
			}
			output++;
			if (output%10==0) newf << endl << "\t\t\t";
		}
	}
	newf << ";\n\nEnd;\n";
	newf.close();
	
	cout << endl;
	cout << "CHARSET for gappy sites has been generated.\n\n";
	cout << endl << "Press return to quit..." << endl;
	cin.get();

	return 0;
	}
