// SPECIES RICHNESS ANALYSIS
// Program will read from text file species and their abundance at each site.
// The expected random distribution of species will then be calculated based on 10000 iterations.
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <limits.h>
#include <ctime>
#include <string>
#include <math.h>
#include <sys/timeb.h>


#define SPP 999			//Maximum number of species = 999
#define NAMLEN 30		// Species identifiers can have no more than 30 characters
#define SITE 200		//The  maximum number of samples allowed.
						// NOTE that throughout code I use "site" as equivalent to "sample".
#define ITER 10000	// number iterations to run
#define SITEGP 22
double ran2(long * idum);


int main()  {

using namespace std;

long spp;
int j=0;
int m=0;
int lastm=0;
int startsite;
long number=0;
char name[SPP][NAMLEN];		//variable that holds species names
long abund[SPP];			//variable that holds abundance values for each species (read from file)
double array[SPP][SITE];	//array to keep track of the abundance of each species in each sample
long fill[SITE];			//variable that holds abundance values for each sample (user input)
long group;
long groups =0;
long sitenumber;
long sitegroup[SITEGP];		//variable to keep track of which group sample belongs
long temp2;
long actsites=0;
double rannum=0;
long i, k, z, num, site, ind, count;
double ranind=0;
long bee=0;
long iter;
long actspp = 0;
long totind = 0;
long totind2 = 0;
int sitecheck = 0;
char response;
string filein;
string fileout;
char buffer[32];
string label0;
string label1;
string label2;
int outbuf;
int pos;


long *idum;
long first;
struct timeb timebuffer;		//used for random number generation

ftime(&timebuffer);
idum=&first;
first = -timebuffer.time + (timebuffer.millitm*100000);



/*open file of data*/
fprintf(stdout, "What is the name of the file with data? (no more than 20 characters)\n");
fprintf(stdout, "File should be plain text in the format of:\n");
fprintf(stdout, "\nspecies 1 \nabundance of species 1 \nspecies 2\nabundance of species 2\netc.");
fprintf(stdout, "\n\n for all species in your data set\n");

cin >> filein;
//fprintf(stdout, "%s \n", filein);

fprintf(stdout, "Enter a file name for output (no more than 20 characters).\nIF FILE EXISTS, IT WILL BE ERASED.\n");
cin >> fileout;
cin.get();

//read in data from file

ifstream fin;
fin.clear();
fin.open(filein.c_str());
	if(!fin.is_open()) {
		cerr << "ERROR: Data input file could not be openned.  Program terminated.  Hit enter.\n";
		cin.get();
                exit(1);
		}


//read in data
if (fin.bad()) {cerr << "error in file stream. hit enter.\n"; cin.get(); exit(1);}

//std::getline(fin, label0, '\n');

fin >> label0;
if(label0.compare("samples")==0){
            while (!(fin >> actsites) || (actsites>700)) {
                cerr << "ERROR in reading data for number of samples. Program terminated.  Hit enter.\n";
                cin.get();
                fin.clear();
                fin.close();
                exit(1);
                }
	for(k=0; k<actsites; k++) {
            fin >> fill[k];
		if ( (fin.rdstate() & ifstream::failbit ) != 0 ) {
			cerr << "ERROR in reading data for sample abundance " << k+1 <<". Program terminated.  Hit enter.\n";
			cin.get();
			 fin.clear();
			 fin.close();
			 exit(1);
                }
  //          cout << fin.tellg() << endl;
//            fin.clear();
            totind2+=fill[k];
        }
}




pos=fin.tellg();
fin >> label1;
if(label1.compare("groups")==0) {
	while(!(fin >> groups) || (groups>actsites)) {
		cerr << "ERROR in number of groups read from file.\nMust be a number less than the number of samples.\nProgram terminated.  Hit enter.\n";
		cin.get();
		fin.clear();
		fin.close();
		exit(1);
	}
	for(k=0; k<groups; k++) {
		while(!(fin >> sitegroup[k])) {
			cerr << "ERROR in reading number of samples in group " << (k+1) <<".\nProgram terminated. Hit enter.\n";
			cin.get();
			fin.clear();
			fin.close();
			exit(1);
		}
//		fin.ignore();
	}
}
else fin.seekg(pos);


fin >> label2;
if(label2.compare("species")==0) {
		for(spp=0;spp<SPP && !fin.eof();spp++) {
		  fin >> name[spp];
		  fin >> abund[spp];
		  actspp++;			//count number of species read in
		}
}
else {
	cerr << "No 'species' label in input file.  Program terminated.  Hit enter. \n";

        cin.get();
	exit(1);
}
fin.clear();
fin.close();


// Count up number of total # of individuals read in
for (i=0; i<actspp; i++) {
	totind += abund[i];
}

cout << actspp << "  species read from file." << endl;
cout << "total abundance for all species combined is " << totind << endl;



//CHECK FOR BAD DATA ENTRY
while (totind2 != totind) {
	cout << "Error - your sample abundances do not match the overall abundance" <<endl;
	cout << "Sum of sample abundances = " << totind2 << endl;
	cout << "Sum of species abundances = " << totind << endl;
	cout << "\nPlease correct file data and re-run program" << endl;
	cout << "Hit enter to quit program."<< endl;
	cin.get();
	exit(1);

}




//WRITE COLUMN HEADERS to output file

	// open file for writing
	fstream fout;
	fout.open(fileout.c_str(), ios::out | ios_base::trunc);
			if(!fout.is_open()) {
			cerr << "couldn't open output file\n";
			}

	for(site=0;site<actsites;site++)  {
    	fout << "SAMPLE_"  << site+1 << " ";
		}
	for(group=0; group<groups;group++) {
		fout << "GROUP_" << group+1 << " ";
		}
	fout << "\n";
	fout.clear();
	fout.close();


// START ITERATIONS

for(iter=0;iter<ITER;iter++)  {
	cout << iter << " ";

// Randomize individuals among sites

// first copy data into temp array vector
//temp has totind items, each item is a spp number, indicating which species that ind. is
z=0;
long *temp = new long [totind];

	for(spp=0;spp<actspp;spp++) {
		num=abund[spp];
		while (num>0) {
		  temp[z] = spp;
		  num--;
		  z++;
		}
	}

// 12 Dec 04 - checked printout of above and data transferred correctly



//Create pseudo_data arrays and set to all zeroes

long *spprich = new long [actsites];
long *sitecnt = new long [actsites];
long *gpspprich = new long [(groups+3)]; //add 3 b/c if no groups, value =-2 & error caused

		for(spp=0; spp<actspp; spp++) {
			for(site=0; site<actsites; site++) {
				array[spp][site]=0;
			}
		}
		for(site=0;site<actsites;site++) {
			spprich[site]=0;  }


		for(k=0; k<groups; k++) {
			gpspprich[k]=0;  }


//set counting parameters to zero
	for(site=0;site<actsites;site++) {
		sitecnt[site]=0;
	}



count=totind;
// now assign individuals at random to array

		while (count > 0) {			//continue until all individuals assigned
		ranind = ran2(idum)*totind;		//choose random individual
		bee = (int) ranind;
		rannum = ran2(idum)*actsites;	//choose random sample for individual to go to
		number = (int) rannum;
		spp = temp[bee];
		site = number;
		if (sitecnt[site] < fill[site] && spp != -99) {  //check that individual has not already been placed and that sample choosen not already full
			array[spp][site]++;
			temp[bee]=-99;
			}
		else continue;
	  count--;
	  sitecnt[site]++;
	}




//Re-set arrays of species richness to zero in prep. for this iteration's calculation
//for (j=0; j<actsites;j++) {
//	spprich[j]=0; }

/*//DEBUG PRINT OF ARRAY DATA
	ofstream fout2;
	fout2.open(fileout, ios::app);
			if(!fout2.is_open()) {
			cerr << "couldn't open output file\n";
			}

		for (site=0; site<actsites; site++) {
			for(spp=0; spp<actspp;spp++) {
			fout2 << site << " " << spp << " " << array[spp][site] << endl;
			}
		}
	fout2 << "\n";
	fout2.clear();
	fout2.close();
*/

//Summary stats for each site separately
	for (site=0; site<actsites; site++) {
		for(spp=0; spp<actspp;spp++) {
			if(array[spp][site] > 0) spprich[site]++;
		}
	}
//Summary stats for site groups
	startsite=0;

	if(groups>0) {
		for(k=0;k<groups;k++) {
				sitenumber=sitegroup[k];
				temp2=0;
				for (spp=0;spp<actspp;spp++) {
					temp2=0;
					for(m=startsite; m<(startsite+sitenumber); m++) {
						if(array[spp][m]>0) temp2++;
					}
				if(temp2>0) gpspprich[k]++;
				}
		  startsite+=sitenumber;
		}

	}

delete [] temp;  //free memory for dynamic arrays
delete [] sitecnt;


// ouput randomized data to file

// open file for writing
	ofstream fout;
	fout.open(fileout.c_str(), ios::out | ios_base::app);
			if(!fout.is_open()) {
			cerr << "couldn't open output file\n";
			exit(1);
			}

	for(site=0;site<actsites;site++)  {
    	fout << spprich[site] << " ";
		}
	for(group=0; group<groups;group++) {
		fout << gpspprich[group] << " ";
	}
	fout << "\n";
	fout.clear();
	fout.close();

delete [] spprich;
delete [] gpspprich;
} // end interations loop



return 0;
}