/*----------------------------------------------------------------------------------------------------------------------
| The "reorder" program by Edward L. Braun, Dept. of Biology, Univ. Florida.
| Modified version of Paul Lewis' ncltest program that will rearrange taxa in user specified orders
|   distribution completed in 2009
*/
#include <string.h>
#include "ncl.h"
#include "ncltest.h"

#define MAXSITE 100000
#define MAXTAXA 5000
#define MAXTAXSETINFO 100000

/*----------------------------------------------------------------------------------------------------------------------
|	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[])
	{
	char infname[256];
	char outfname[256];
	char newfname[256];
	char authfname[256];
	char missingfname[256];
	char scratchfname[256];
	char charsetfname[256];
	char taxsetfname[256];
	char finalcharsetfname[256];
	char finaltaxsetfname[256];
	bool outputnexus=true;
	unsigned offset=0; // one less than the first character of the alignment
	strcpy(infname, "infile.nex");
	strcpy(outfname, "output.txt");
	strcpy(newfname, "outfile.nex");
	strcpy(authfname, "authority.txt");
	if (argc==1) {
	  cerr << "Reorder taxa based upon authority file. Usage:\n";
	  cerr << "  $ ./reorder <infile> <authority> <outfile> <offset>\n";
	  cerr << "\n  All command line arguments are optional.\n";
	  cerr << "  this program will default to the following input and output files:\n";
	  cerr << "    infile.nex\n";
	  cerr << "    authority.txt\n";
	  cerr << "    outfile.nex\n";
	  cerr << "    offset - used if called by the supernexus wrapper.\n";
	  cerr << "  This program also creates a file called 'output.txt' that records the\n";
	  cerr << "  contents of the nexus input file.\n";
	}
	if (argc > 1)
		strcpy(infname, argv[1]);
	if (argc > 2)
		strcpy(authfname, argv[2]);
	if (argc > 3)
		strcpy(newfname, argv[3]);
	if (argc > 4) {
		offset=atoi(argv[4])-1;
		outputnexus=false;
	}
	strcpy(missingfname,"missing.");
	strcat(missingfname,newfname);
	strcat(missingfname,".txt");

	/*----------------------------------------------------------------------------------------------------------------------
	|  Read the authority file, including the name of a single set to include.
	*/
	ifstream authf(authfname);
	char tempname[256];
	char authname[MAXTAXA][256];
	bool includeonly=false;
	bool includeonlyfound=false;
	char includeonlysetname[256];
	char includeonlysites[MAXSITE];
	tempname[0]='-';
	unsigned nauthtax=0;
	cerr << "Reading the authority file:\n";
	while (tempname[0]!='\0') {
		tempname[0]='\0';
		authf >> tempname;
		if (tempname[0]=='#' && nauthtax==0) {
			// indicate a charset that will be the only set included
			tempname[0]='\0';
			authf >> tempname;
			if (tempname[0]=='\0') {
				cerr << "Authority file does not include taxon names.\n";
				return 0;
			}
			includeonly=true;
			strcpy(includeonlysetname,tempname);
			cerr << "\nInclude only characters in CHARSET " << includeonlysetname << ".\n\n";
			tempname[0]='\0';
			authf >> tempname;
		}
		if (tempname[0]!='\0') {
			strcpy(authname[nauthtax],tempname);
			nauthtax++;
		}
	}
	if (outputnexus) cerr << "There are a total of " << nauthtax << " taxa in the authority file.\n";
	authf.close();

	char setinfo[MAXSITE]; // String used to hold set information

	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);
	
	// Variables used for analyses
	unsigned i, j, k, states, fflag, eflag;
	unsigned width = taxa->GetMaxTaxonLabelLength();
	int ntax = characters->GetNTax();
	cerr << "Ntax  = " << ntax << endl;
	int nchar = characters->GetNChar();
	cerr << "Nchar = " << nchar << endl;

	NxsStringVector z;
	
	unsigned ntaxsets = assumptions->GetNumTaxSets();
	strcpy(taxsetfname,"scratch.rawtaxsetinfo.");
	strcat(taxsetfname,infname);
//	strcpy(charsetfname,".txt");
	ofstream taxscrf(taxsetfname);
	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) << ": ";
			NxsString temptaxsetname(*zi);
			taxscrf << temptaxsetname.BlanksToUnderscores() << " ";
			NxsUnsignedSet &x = assumptions->GetTaxSet(*zi);
			NxsUnsignedSet::const_iterator xi;
			for (xi = x.begin(); xi != x.end(); xi++) {
				nexus.outf << (*xi + 1) << " ";
				taxscrf << characters->GetTaxonLabel(*xi) << " ";
			}
			nexus.outf << endl;
			taxscrf << ";\n";
			}
		}
	taxscrf.close();

	unsigned ncharsets = 0;
	strcpy(charsetfname,"scratch.rawcharsetinfo.");
	strcat(charsetfname,infname);
//	strcpy(charsetfname,".txt");
	ofstream charscrf(charsetfname);

	NxsString includeonlyNexsetname(includeonlysetname); // used if "includeonly" is true
	for (i=0; i<nchar; i++) {
		includeonlysites[i]='1';
		includeonlysites[i+1]='\0';
	}
	if (assumptions->GetNumCharSets() > 0)
		{
		assumptions->GetCharSetNames(z);
		NxsStringVector::const_iterator zi;
		for (zi = z.begin(); zi != z.end(); zi++)
			{
			NxsString tempsetname(*zi);
			nexus.outf << "Charset " << tempsetname.BlanksToUnderscores() << ": ";
		//	nexus.outf << "Charset " << (*zi) << ": ";
			charscrf << tempsetname.BlanksToUnderscores() << " ";
			// Initialize the string used to hold set information
			for (i = 0; i < nchar; i++) setinfo[i]='0';
			setinfo[nchar]='\0';
			NxsUnsignedSet &x = assumptions->GetCharSet(*zi);
			NxsUnsignedSet::const_iterator xi;
			for (xi = x.begin(); xi != x.end(); xi++) {
				nexus.outf << (*xi + 1) << " ";
				setinfo[(*xi)]='1';
				if (includeonly && tempsetname.Equals(includeonlyNexsetname)) {
					includeonlysites[(*xi)]='0';
				}
			}
			nexus.outf << endl;
			charscrf << setinfo << endl;
			ncharsets++;
			}
		}
//	cerr << includeonlysites << endl;
		
	charscrf.close();

	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;
			}
		}

	/*----------------------------------------------------------------------------------------------------------------------
	|  Output a nexus file with the same order as the authority file.
	*/

	strcpy(scratchfname,"scratch.");
	strcat(scratchfname,infname);
	ofstream scratchf(scratchfname);
	/*----------------------------------------------------------------------------------------------------------------------
	|  Identify the sites that are all gaps or missing.
	*/
	char allmissing[MAXSITE];
	for (i = 0; i < nchar; i++) allmissing[i]='1';
	allmissing[nchar]='\0'; // Set up a string with information about whether sites are all gaps or missing
	/*----------------------------------------------------------------------------------------------------------------------
	|  Output the rearranged data block.
	*/
//	scratchf << "#NEXUS\n\nBegin data;\n\tDimensions ntax=";
//	scratchf << ntax << " nchar=" << nchar << ";" << endl;
//	scratchf << "\tFormat datatype=dna gap=- missing=?;\n\tmatrix\n";
	cerr << "Output order:\n";
	char datum; char ambigdatum[256]; eflag=0;
	ofstream missingf(missingfname);                // missingf will have include a taxset with those taxa in the authority
	missingf << "\tTaxset No_" << infname << " = "; // file that were not found in the input file.
	for (j = 0; j < nauthtax; j++)
		{
		fflag=0;
	//	authf >> authname[j];
		cerr << "   " << authname[j];
		NxsString authTaxonLabel(authname[j]);
		for (i = 0; i < ntax; i++)
			{
			//  Read the currentTaxonLabel into a NxsString...
			NxsString currTaxonLabel = characters->GetTaxonLabel(i);
			// determine whether the currTaxonLabel matches the TaxonLabel from the authority file
			currTaxonLabel.BlanksToUnderscores();
			if (currTaxonLabel.Equals(authTaxonLabel)) {
				fflag=1; // set fflag (found flag) to 1
				cerr << " --> taxon found" << endl;
			
				scratchf << currTaxonLabel;
				scratchf << " ";

				// Print out enough spaces to even up the left edge of the matrix output
				//
			//	unsigned currTaxonLabelLen = currTaxonLabel.size();
			//	unsigned diff = width - currTaxonLabelLen;
			//	for (k = 0; k < diff+5; k++) {
			//		scratchf << ' ';
			//	}  // end for (k = 0...
	
				for (unsigned currChar = 0; currChar < nchar; currChar++)
					{
					if (characters->IsGapState(i, currChar)) {
						datum='-';
						scratchf << datum;
					} // end if (characters->IsGapState(i, currChar)...
					else if (characters->IsMissingState(i, currChar)) {
							datum='?';
							scratchf << datum;
						} // end if (characters->IsMissingState(i, currChar)...
						else {
							allmissing[currChar]='0'; // '0' indicates that at least one taxon is defined (not gap or missing) 
							states=characters->GetNumStates(i, currChar);
							if (states==1) {
								scratchf << characters->GetState(i, currChar, 0);
							}  // end if (states==1...
							else {
							for (int z=0; z<states; z++) 
								{
									ambigdatum[z]=characters->GetState(i, currChar, z);
									ambigdatum[z+1]='\0';
								} // end for (int z=0...
								scratchf << ambigDNA(ambigdatum,states);
							} // end else
						} // end else
					} // for (unsigned currChar = 0...
					scratchf << endl;					
			} // if (currTaxonLabel.Equals(authTaxonLabel)...
		} // for (i = 0; i < ntax; i++...
		if (includeonly) {
			for (i = 0; i < nchar; i++) {
				if (includeonlysites[i]=='1') allmissing[i]='1';
			}
		} // end if (includeonly...
		if (fflag==0) {
			cerr << " --> ** ERROR ** Taxon was not found" << endl;
			scratchf << authTaxonLabel;
			scratchf << " ";
			for (unsigned dummyChar = 0; dummyChar < nchar; dummyChar++) scratchf << "?";
			scratchf << endl;
			missingf << authTaxonLabel << ' ';
			eflag++;
		} // if (fflag==0...
	} // for (j = 0; j < ntax; j++...
//	scratchf << "\t;\nEnd;\n\n";
	scratchf.close();
	missingf << " ;\n";
	missingf.close();
	
	cerr << endl;
	cerr << "File has been rearranged.";
	if (eflag>0) {
		cerr << " " << eflag << " taxa in the authority file were not found.";
	}
	cerr << endl;
	
	/*----------------------------------------------------------------------------------------------------------------------
	|  Determine the number of sites that were retained in the alignment (after eliminating columns that are all gaps or
	|  missing given the taxa included). Use this to modify the set information.
	*/
	unsigned nincluchar=0;
	for (k = 0; k < nchar; k++) {
		if (allmissing[k]=='0') nincluchar++;
	}
	char setname[256];
	char ch;
	strcpy(finalcharsetfname,"scratch.charsetinfo.");
	strcat(finalcharsetfname,infname);
//	strcpy(finalcharsetfname,".txt");
	ifstream rawtaxscrf(charsetfname);
	ofstream outtaxscrf(finalcharsetfname);
	outtaxscrf << ncharsets << endl; // first output the number of charsets
	for (i = 0; i < ncharsets; i++) {
		rawtaxscrf >> setname;
		outtaxscrf << setname << " ";
		for (k = 0; k < nchar; k++) {
			rawtaxscrf >> ch;
			if (allmissing[k]=='0') outtaxscrf << ch;
		}
		outtaxscrf << endl;
	}
	rawtaxscrf.close();
	outtaxscrf.close();

	/*----------------------------------------------------------------------------------------------------------------------
	|  Either output the header to a nexus file and the rearranged sequence data as a nexus data block or simply
	|  output the rearranged sequence data.
	*/
	char taxname[256];
	char nt;
	unsigned currTaxonLen, diff;
	ifstream iscrf(scratchfname);
	ofstream newf(newfname);
	if (outputnexus) newf << "#NEXUS\n\nBegin data;\n\tDimensions ntax=";
	if (outputnexus) newf << nauthtax << " nchar=";
	if (!outputnexus) newf << offset+1 << '\t';
	newf << nincluchar+offset;
	if (outputnexus) newf << ";";
	newf << endl;
	if (outputnexus) newf << "\tFormat datatype=dna gap=- missing=?;\n\tmatrix\n";
	for (j = 0; j < nauthtax; j++) {
		iscrf >> taxname;
		newf << taxname;
		// Print out enough spaces to even up the left edge of the matrix output
		currTaxonLen = strlen(taxname);
		diff = width - currTaxonLen;
		for (k = 0; k < diff+5; k++) {
			newf << ' ';
		}  // end for (k = 0...
		// Now output the matrix
		for (k = 0; k < nchar; k++) {
			iscrf >> nt;
			if (allmissing[k]=='0') {
				newf << nt;
			}
		}  // end for (k = 0...
		newf << endl;
	}
	if (outputnexus) newf << "\t;\nEnd;";
	newf << "\n\n";
	iscrf.close();
	newf.close();
	
	/*----------------------------------------------------------------------------------------------------------------------
	|  Output the number of charsets to a second scratchfile if 'outputnexus' is false.
	*/
	if (!outputnexus) {
		strcat(newfname,".sets");
		ofstream newheaderf(newfname);
		newheaderf << ncharsets << endl;
		newheaderf.close();
	}
	
	/*----------------------------------------------------------------------------------------------------------------------
	|  Output the charset and taxset information as a nexus assumptions block.
	*/
	unsigned posit, p, printed, plen;
	bool isemptyset;
	ofstream newsetf(newfname,ios::app);
	ifstream incharscrf(finalcharsetfname);
	incharscrf >> ncharsets;
	if (outputnexus) newsetf << "\nBegin assumptions;\n";
	for (i = 0; i < ncharsets; i++) {
		incharscrf >> setname;
		cerr << "Checking CHARSET " << setname << "...";
		incharscrf.getline(setinfo,MAXSITE,'\n');
		plen=strlen(setinfo);
		// check whether the set is an empty set
		p=0;	isemptyset=true;
		while (isemptyset && p<plen) {
			if (setinfo[p] == '1') isemptyset=false;
			p++;
		} // end while (p<plen...
		if (!isemptyset) {
			cerr << "not empty\n";
			newsetf << "\tCharset ";
			newsetf << setname << " = ";
			p=0;	printed=0;
			while ( p<plen ) {
				if ( setinfo[p] == '1' ) {
					if ( setinfo[p+1] == '1' && setinfo[p+2] == '1' ) {
						newsetf << p+offset << " - ";
						while ( setinfo[p]=='1' ) { p++; }
						newsetf << p+offset-1 << ' ';
						printed++;
					} // end if ( setinfo[p+1] == '1' && setinfo[p+2] == '1' ...
					else {
						newsetf << p+offset << ' ';
						printed++;
					} // end else
					if (printed==20) {
						newsetf << endl << "\t\t\t";
						printed=0;
					} // end if (printed==20...
				} // if ( setinfo[p] == '1'...
				p++;
			} // while ( p<plen...
			newsetf << ";\n";
		} // end if (!isemptyset...
		else cerr << "EMPTY\n";
	} // end for (i = 0; i < ncharsets...
	incharscrf.close();
	
	// finaltaxsetname is used here to facilitate further modification
	strcpy(finaltaxsetfname,"scratch.rawtaxsetinfo.");
	strcat(finaltaxsetfname,infname);
//	strcpy(finaltaxsetfname,".txt");
	ifstream intaxscrf(finaltaxsetfname);
	char taxsetinfo[MAXTAXSETINFO];
	bool emptytaxset;
	for (i = 0; i < ntaxsets; i++) {
		emptytaxset=true;
		intaxscrf >> setname;
		strcpy(taxsetinfo,"Taxset ");
		strcat(taxsetinfo,setname);
		strcat(taxsetinfo," = ");
	//	newsetf << "\tTaxset " << setname << " = ";
	//	printed=0;
		while (setname[0]!=';') {
			intaxscrf >> setname;
			for (j = 0; j < nauthtax; j++) {
				if (setname[0]==';') {
					strcat(taxsetinfo,";");
					j=nauthtax;
				}
				if (strcmp(setname,authname[j])==0) {
					emptytaxset=false;
					strcat(taxsetinfo,setname);
					strcat(taxsetinfo," ");
					j=nauthtax;
				}
			}
		//	newsetf << setname << " ";
		//	if (printed==10) {
		//		newsetf << endl << "\t\t\t";
		//		printed=0;
		//	} // end if (printed==20...
		} // end while (setname[0]!=';'...
		if (!emptytaxset) newsetf << '\t' << taxsetinfo << endl;
	} // end for (i = 0; i < ntaxsets...
	if (outputnexus) newsetf << "End;\n";
	intaxscrf.close();
	newsetf.close();
	
	if (outputnexus) {
		system("rm scratch.*");
		cerr << endl << "Press return to quit..." << endl;
		cin.get();
	}

	return 0;
}
