/*
 
 Copyright (C) 2010,2011 Timo Lassmann <timolassmann@gmail.com>
 
 This file is part of aln_filter.
 
 aln_filter is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 aln_filter is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with aln_filter.  If not, see <http://www.gnu.org/licenses/>.
 
 */
#include "aln_filter.h"
#include "interface.h"
#include <stdio.h>

int main (int argc,char * argv[])
{
	struct db* db = 0;
	
	struct parameters* param = 0;
	param =  interface( param,argc, argv);
	// insert code here...
	//printf("Hello, World!\n");
	if(!param->outfile){
		outfile = stdout;
	}else{
		if ((outfile = fopen(param->outfile, "w")) == NULL){
			fprintf(stderr,"can't open output\n");
			exit(-1);
		}
	}
	db = fix_sam_file(db,param);
	if(param->outfile){
		fclose(outfile);
	}
	
	free_param(param);
	free_db(db);
	return 0;
}

void free_db(struct db* db)
{
	int i; 
	for(i = 0; i < db->seq_info->num_seq;i++){
		free(db->seq_info->names[i]); 
	}
	free(db->dbname);
	free(db->seq_info->names);
	free(db->seq_info->sn_len);
	free(db->seq_info->len);
	free(db->seq_info->nuc_counts);
	free(db->seq_info);
	free(db->bfa);
	free(db);
}

struct db* fix_sam_file(struct db* db,struct parameters* param)
{
	FILE* file = 0;
	char tmp_text[MAX_LINE];
	char line[MAX_LINE];
	char command[1000];
	int i,j;
	int tmp = 0;
	unsigned int pos;
	unsigned int aln_len;
	unsigned int add = 0;
	unsigned char* real_chr_name;
	unsigned int chr_len = 0;
	
	char readname[100];
	int map_flag;
	char chromosome[100];
	int read_pos;
	//int aln_len = 0;
	int quality;
	char cigar[100];
	char rnext[100];
	char pnext[100];
	char TLEN[100]; 
	char seq[1000];
	char base_quality[1000];
	char remainder[1000];
	int errors; 

	int storage_num_chr = 10; 
	db = malloc(sizeof(struct db));
	
	db->seq_info = malloc(sizeof(struct db_seq_info));
	db->seq_info->names = malloc(sizeof(char*)* storage_num_chr);
	db->seq_info->sn_len = malloc(sizeof(int)* storage_num_chr);
	db->seq_info->len = malloc(sizeof(unsigned int)* storage_num_chr);
	db->seq_info->nuc_counts = malloc(sizeof(float)*5);
	
	db->bfa = 0;
	db->dbname  = malloc(sizeof(char)*(2));
	db->dbname[0] = 'A';
	db->dbname[1] = 0;
	
	int read_seq_names = 0;
	
	unsigned int real_start = 0;
	
	int sam = 0;

	if(!strcmp(".sam", param->infile[0] + (strlen(param->infile[0] ) - 4))){
		sam = 1;
		//}else if (byg_end(".bam", param->infile[i])  == strlen(param->infile[i])){
	}else if (!strcmp(".bam", param->infile[0] + (strlen(param->infile[0] ) - 4))){
		sam = 2;
		//}else if (byg_end(".fa", param->infile[i])  == strlen(param->infile[i])){
	}
	
	if(!sam){
		fprintf(stderr,"only works with sam / bam files.%s \n", param->infile[0]);
		exit(-1);
	}
	if(sam == 2){
		command[0] = 0;
		strcat ( command, "samtools view -H ");     // Copy name into full name
		i = sprintf (tmp_text, "%s ", param->infile[0]);
		strcat ( command, tmp_text);
		//fprintf(stderr,"%s\n",command);
		if (!(file = popen(command, "r"))) {
			fprintf(stderr,"Cannot open bam file '%s' with command:%s\n",param->infile[0],command);
			exit(-1);
		}
	}else if(sam == 1){
		command[0] = 0;
		strcat ( command, "samtools view -SH ");     // Copy name into full name
		i = sprintf (tmp_text, "%s ", param->infile[0]);
		strcat ( command, tmp_text);
		//fprintf(stderr,"%s\n",command);
		if (!(file = popen(command, "r"))) {
			fprintf(stderr,"Cannot open sam file '%s' with command:%s\n",param->infile[0],command);
			exit(-1);
		}
	}else{
		if (!(file = fopen(param->infile[0] , "r" ))){
			fprintf(stderr,"Cannot open sam file '%s'\n",param->infile[0]);
			exit(-1);
		}
		
	}
	
	while(fgets(line, MAX_LINE, file)){
		fprintf(outfile,"%s",line);
		if(line[0] == '@'){
			if(line[1] == 'S' && line[2] == 'Q' ){
				
				tmp = byg_end("SN:", line );
				if(tmp){
					db->seq_info->names[read_seq_names] = malloc(sizeof(unsigned char)*(50) );
					add = 0;
					for(j = tmp;j < MAX_LINE;j++){
						if(isspace((int)line[j])){
							db->seq_info->names[read_seq_names][add] = 0;
							break;
						}
						db->seq_info->names[read_seq_names][add] = line[j];
						add++;
					}
				}
				
				//db->seq_info->names[read] = (unsigned char*)strdup(name);
				db->seq_info->sn_len[read_seq_names] = add;
				tmp = byg_end("LN:", line );
				db->seq_info->len[read_seq_names] = atoi(line+tmp);
				read_seq_names++;
				if(read_seq_names == storage_num_chr){
					storage_num_chr = storage_num_chr << 1;
					db->seq_info->names = realloc(db->seq_info->names,sizeof(char*)* storage_num_chr);
					db->seq_info->sn_len = realloc(db->seq_info->sn_len,sizeof(int)* storage_num_chr);
					db->seq_info->len = realloc(db->seq_info->len,sizeof(int)* storage_num_chr);
				}
			}
			
		}else{
			break;
		}
	}
	if(!read_seq_names){
		fprintf(stderr,"SAM header seems to be absent.\n");
		exit(-1);
		fclose(file);
		return db;
		
	}
	
	db->seq_info->num_seq = read_seq_names;
	
	
	
	//real_start = 0;
	pos = 0;
	for(i = 0; i < db->seq_info->num_seq ;i++){
		real_start = db->seq_info->len[i];
		db->seq_info->len[i] = pos;
		pos = real_start + pos;
	}
	//real_start = 0;
	int equal;
	if(sam == 2){
		pclose(file);
		command[0] = 0;
		strcat ( command, "samtools view  ");     // Copy name into full name
		i = sprintf (tmp_text, "%s ", param->infile[0]);
		strcat ( command, tmp_text);
		//fprintf(stderr,"%s\n",command);
		if (!(file = popen(command, "r"))) {
			fprintf(stderr,"Cannot open bam file '%s' with command:%s\n",param->infile[1],command);
			exit(-1);
		}
	}else if(sam == 1){
		pclose(file);
		command[0] = 0;
		strcat ( command, "samtools view -S ");     // Copy name into full name
		i = sprintf (tmp_text, "%s ", param->infile[0]);
		strcat ( command, tmp_text);
		//fprintf(stderr,"%s\n",command);
		if (!(file = popen(command, "r"))) {
			fprintf(stderr,"Cannot open sam file '%s' with command:%s\n",param->infile[1],command);
			exit(-1);
		}
	}else{
		rewind(file);		
	}
	
	
	while(fgets(line, MAX_LINE, file)){
		//VHE-229301321014-18-841-1-9426	0	chr1	10070	1	8M1I33M	*	0	0	CCCTAACCTCTAACCCTAACCCTAACCCTAACCCTAACCCAA	*	NM:i:1	MD:Z:41	XP:Z:~33:::2~!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		
		sscanf(line, "%s	%d	%s	%d	%d	%s	%s	%s	%s	%s	%s	%[^\n]",readname,&map_flag, chromosome,&read_pos,&quality, cigar ,rnext,pnext ,TLEN,seq,base_quality,remainder);
		
		aln_len =get_aln_length(cigar); 
		//fprintf(stderr,"%s\n%d\n",line,get_aln_length(cigar));
		errors = 0;
		tmp = byg_end("NM:i:", remainder  );
		if(tmp){
			errors = (float) atoi(remainder+tmp);
		}else{
			errors= -1;								
		}	
		equal = 0;
		for(chr_len = 0;chr_len < strlen(seq);chr_len++){
			if(seq[chr_len] == '='){
				equal++;
			}
		}
		
		if(chromosome[0] != '*'){
			
			//if(aln_len - errors != equal){
			
			//	fprintf(stderr,"%d	%d	%d	%f\n",errors,aln_len, equal, ((float) aln_len  -  (float)(errors)) / (float) aln_len);
			//	fprintf(stderr,"%f	%s",1.0  - ((float) aln_len  -  (float)(errors)) / (float) aln_len,  line);
			//}
			if(errors == -1){
				fprintf(stderr,"WARNING: SAM line does not contain the NM:i: field.\n");
				fprintf(outfile,"%s",line);
			}else {
					
				if(((float) aln_len  -  (float)(errors)) / (float) aln_len < param->id_cutoff){
					map_flag = map_flag | 0x200;
				}
				
				chr_len = get_chr_offset(db->seq_info, chromosome);
				if(chr_len == 0xFFFFFFFFu){
					fprintf(stderr,"ERROR: %s not found in indexed sequences\n",chromosome);
					fprintf(stderr,"%s\n",line);
					exit(-1);
				}
				
				
				//errors =  (float) ri->len / ( 1.0+  ri->identity[0]) ;
				//break;
				//fprintf(stderr,"ERRORS:%d\n\n",errors);
				//pos = chr_len + pos;
				
				if(strcmp((char*)get_chr_name( db->seq_info ,chr_len+read_pos ) , (char*)get_chr_name( db->seq_info ,chr_len+read_pos+aln_len ))){
					fprintf(stderr,"error: alignment crossed chromosome boundaries\n%s",line);
					map_flag = map_flag | 0x200;
					if(!param->discard){
						fprintf(outfile,"%s\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",readname,map_flag, chromosome,read_pos,quality, cigar ,rnext,pnext ,TLEN,seq,base_quality,remainder);
					}
				}else if(strcmp((char*)get_chr_name( db->seq_info ,chr_len+read_pos ) , chromosome)){
					//unsigned char* get_chr_name(struct db_seq_info* si, unsigned int pos)
					real_chr_name = get_chr_name(db->seq_info,chr_len+read_pos);
					fprintf(stderr,"error: should be %s	%d\n%s",real_chr_name ,get_chr_pos(db->seq_info, chr_len+read_pos)+1,line);
					if(!param->discard){
						fprintf(outfile,"%s\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",readname,map_flag, real_chr_name,get_chr_pos(db->seq_info,chr_len+read_pos)+1,quality, cigar ,rnext,pnext ,TLEN,seq,base_quality,remainder);
					}
					
					//fprintf(stderr,"%s	%s\n",base_quality,remainder);
				}else{
					fprintf(outfile,"%s\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",readname,map_flag, chromosome,read_pos,quality, cigar ,rnext,pnext ,TLEN,seq,base_quality,remainder);
				}
			}
		}else{
			fprintf(outfile,"%s",line);
		}
	}
	if(sam==2){
		pclose(file);
	}else if(sam == 1){
		pclose(file);
	}else{
		fclose(file);
	}
	return db;
}



int byg_end(const char* pattern,const char*text)
{
	const char* tmp = 0;
	int Tc;
	int i  = 0;
	int s = 0;
	int T[256];
	for (i = 0;i < 256;i++){ 
		T[i] = 0; 
	}
	
	int m = (int)strlen(pattern);
	int n = (int)strlen(text);
	if (m > n){
		i = m;
		m = n;
		n = i;
		tmp = text;
		text = pattern;
		pattern = tmp;
	}
	
	int mb = (1 << (m-1));
	
	for (i= 0;i < m;i++){
		T[(int)pattern[i]] |= (1 << i);
	}
	
	for (i = 0;i < n;i++){
		s <<= 1;
		s |= 1;
		if(!text[i]){
			return -1;
		}
		Tc = T[(int)text[i]];
		s &= Tc;
		if(s & mb){
			return i+1;
		}
	}
	return 0;
}

int get_aln_length(char* cigar)
{
	char tmp_num[8];
	int l,i,j,c,rp,exit_loop; 
	l = (int) strlen((char*)cigar);
	exit_loop = 0;
	i =0;
	rp = 0;
	//sp = 0;
	while(!exit_loop){
		c = 0;
		if(isdigit((int)cigar[i])){
			j = 0;
			while (isdigit(cigar[i])) {
				tmp_num[j] = cigar[i];
				j++;
				i++;
				if(i == l){
					exit_loop =1;
					break;
				}
			}
			tmp_num[j] = 0;
			
			c = atoi(tmp_num);
		}
		if(isalpha((int)cigar[i])){
			switch (cigar[i]) {
				case 'M':
				case 'I':
				case 'D':
					rp += c;
					break;
				default:
					break;
			}
			c = 0;
		}
		i++;
		if(i == l){
			exit_loop =1;
			break;
		}
		
	}
	
	return rp;
}


unsigned int get_chr_offset(struct db_seq_info* si, char* chr)
{
	int i;
	for(i = 0; i < si->num_seq;i++){
		//fprintf(stderr,"%s	%s\n",chr,si->names[i]);
		if(!strcmp(chr,(char*)si->names[i])){
			return si->len[i];
		}
	}
	return 0xFFFFFFFFu;
}

unsigned char* get_chr_name(struct db_seq_info* si, unsigned int pos)
{
	int i,j;
	j = -1;
	for(i = 1; i< si->num_seq;i++){	
		if(pos < si->len[i]){
			j = i-1;
			break;
		}
	}
	if(j == -1){
		j = si->num_seq-1;
	}
	return si->names[j];
}

unsigned int  get_chr_pos(struct db_seq_info* si, unsigned int pos)
{
	int i,j;
	j = -1;
	for(i = 1; i< si->num_seq;i++){	
		if(pos < si->len[i]){
			j = i-1;
			break;
		}
	}
	if(j == -1){
		j = si->num_seq-1;
	}
	return pos - si->len[j];
}