/***************************************************************************************
ESPRIT: Estimating Species Richness Using Large Collections of 16S rRNA Shotgun Sequences  
Yijun Sun, Yunpeng Cai, Li Liu, Fahong Yu, Michael L. Farrell, William McKendree, William Farmerie
Nucleic Acids Research 2009

Programed by Yunpeng Cai
Copyright Hold by University of Florida, Gainesville, FL32610

****************************************************************************************/
/************************Extracting Sequences From Clustering Results***************/


#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <unistd.h>

#include "FASTA.h"
#include "util.h"

int clnum[20];

#define MAX_SEQS 1000000

int SeqCount=10*MAX_SEQS;

char **SeqStrs;
char **LabelStrs;
char buf[100];
int groupid[10*MAX_SEQS];
FASTA fasta;
char fcname[200];
char cmd[200];
double level;
char prefix[2000];
double leveldn=0.0;
double levelup=1.0;
bool do_trim=false;
int MaxLen=0;

void WriteGroup(double level,int groupidx,int gpptr)
{
		int i,j,slen,ptr;
		#ifdef WIN
		sprintf(fcname,"%s_%.2f\\group_%d.fas",prefix,level,groupidx);
		#else
		sprintf(fcname,"%s_%.2f/group_%d.fas",prefix,level,groupidx);
		#endif
		fasta.OpenWrite(fcname);
		if (do_trim)
		{
			int *pos=(int *)Malloc(MaxLen*sizeof(int));
			char *buf=(char *)Malloc((MaxLen+1)*sizeof(char));
			for (i=0;i<MaxLen;i++) pos[i]=0; 
			for (i=0;i<gpptr;i++)
			{
				slen=strlen(SeqStrs[groupid[i]]);
				for (j=0;j<slen;j++)
				{
					if (SeqStrs[groupid[i]][j] !='-' && SeqStrs[groupid[i]][j] !='_')
						pos[j]=1;
				}
			}
			for (i=0;i<gpptr;i++)
			{
				slen=strlen(SeqStrs[groupid[i]]);
				ptr=0;
				for (j=0;j<slen;j++)
				{
					if (pos[j]>0)
						buf[ptr++]=SeqStrs[groupid[i]][j];
				}
				buf[ptr]=0;
				fasta.WriteFastaSeq(LabelStrs[groupid[i]],buf);
			}	
			free(buf);
			free(pos);
		}
		else
		{
			for (i=0;i<gpptr;i++)
			{
				fasta.WriteFastaSeq(LabelStrs[groupid[i]],SeqStrs[groupid[i]]);
			}
		}	
	fasta.Close();
}

void Usage(){
printf("parsecluster [-t] <fastafile> <clusterfile> [low] [high]\n");
printf("Generate FASTA files of clustered sequences.\n\n");
printf("\tfastafile:\t original or trimmed FASTA sequences.\n");
printf("\tclusterfile:\t corresponding clustering result (.Cluster).\n");
printf("\t[low, high]:\t range of distance levels to generate sequences.\n");
printf("\t\t -t:\t remove common gaps for each group.\n");
}

void GetPrefix(char *output,char *input)
{
	char *pch,*qch;
	if (strchr(input,'.')==NULL)
	{
		strcpy(output,input);
		return;
	}
	pch=input+strlen(input);
	while (pch >input && *pch !='.') 
	{pch--;}
	for (qch=input;qch <pch;qch++)
		*(output++)=*qch;
	*output=0;		
}

int main(int argc, char **argv)
{
	FILE *fp,*fbuf;
	int i,numsq;
	char *label,*seq;
	int rlines=0;
	int bufptr,gpptr,groupidx;
	char gc;
	char c;
	
	while ((c = getopt(argc, argv, "th")) != -1)
    {
        switch (c)
        {
        case 'h':
        	Usage();
        	exit(0);
        	break;
        case 't':
        	do_trim=true;
        	break;
        default:
        	fprintf(stderr,"Unknown option -%c\n",c);
			Usage();
			exit(0);
			break;
        }
    }
	argc -= optind;
	argv += optind;

	
	if (argc <2) 
	{
		Usage();
		exit(0);
	}
	char *seqfile=argv[0];
	char *clfile=argv[1];
	
	if (argc > 2) leveldn=atof(argv[2]);
	if (argc > 3) levelup=atof(argv[3]);
	
	
	SeqStrs=(char **)Malloc(SeqCount*sizeof(char *));
	LabelStrs=(char **)Malloc(SeqCount*sizeof(char *));

	fasta.OpenRead(seqfile);
	i=0;
	
	while (fasta.GetFastaSeq(&label,&seq,0,0))
	{
		LabelStrs[i]=label;
		SeqStrs[i++]=seq;
		if (strlen(seq) >MaxLen)
			MaxLen=strlen(seq);
	}
	fasta.Close();
	numsq=i;

	
	fp=fopen(clfile,"r");
	if (fp==NULL)
	{
		fprintf(stderr,"Cannot Open Cluster File\n");
		exit(0);
	}

	GetPrefix(prefix,clfile);
	printf("Writing Solutions to %s\n",prefix);
	//strcat(prefix,"groups");
	
	while(1)	
	{
		
		level=-1;
		fscanf(fp,"%lf",&level);
		if (level < 0) break;
		if (level >1)
		{
			fprintf(stderr,"Should not get here.\n");
		}
		
		if (level >=leveldn-1e-4 && level <=levelup+1e-4)
		{
			sprintf(cmd,"mkdir %s_%.2f\n",prefix,level);
			printf(cmd);
			system(cmd);
			printf("Level=%.2f\n",level);
		}
		else
		{
			printf("Skipping Level %.2f\n",level);
			do {(gc=fgetc(fp));} while (gc!='\n');
			continue;
		}
		bufptr=0;
		gpptr=0;
		groupidx=0;
		
		while (1)
		{
			gc=fgetc(fp);
			if (gc=='\n' || gc=='\r')
			{
				if (bufptr>0)
				{
					buf[bufptr]='\0';
					groupid[gpptr++]=atoi(buf);
					bufptr=0;
				}

				if (gpptr >0)
				{
					WriteGroup(level,groupidx,gpptr);
					groupidx++;
				}
				printf("Level %.2f Groups %d\n",level, groupidx);
				break;
			}					
			if (gc=='|')
			{
				if (bufptr>0)
				{
					buf[bufptr]='\0';
					groupid[gpptr++]=atoi(buf);
					bufptr=0;
				}

				if (gpptr >0)
				{
					WriteGroup(level,groupidx,gpptr);
					groupidx++;
				}
				gpptr=0;			
			}
			else if (gc==' ')
			{
				if (bufptr>0)
				{
					buf[bufptr]='\0';
					groupid[gpptr++]=atoi(buf);
					bufptr=0;
				}
			}
			else
			{
				buf[bufptr++]=gc;
			}
		}			
	}
	fclose(fp);
	
	for (i=0;i<numsq;i++)
	{
		free(SeqStrs[i]);
		free(LabelStrs[i]);
	}
	free(SeqStrs);
	free(LabelStrs);
	return 0;
}
