// -*- mode: c++; indent-tabs-mode: nil; -*-
//
// Copyright 2009 Illumina, Inc.
//
// This software is covered by the "Illumina Genome Analyzer Software
// License Agreement" and the "Illumina Source Code License Agreement",
// and certain third party copyright/licenses, and any user of this
// source file is bound by the terms therein (see accompanying files
// Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
// Illumina_Source_Code_License_Agreement.pdf and third party
// copyright/license notices).
//
//

//
// makeGStudioBAMIndex.cpp : Creating the Genome Studio BAM linear index
//

/*
 * Author : Mohammed D. Alam
 * This source file implements BAM linear index creation for use in
 * Genome Studio. This implementation depends on the original SAMTools
 * API published at samtools.soureforge.net, and zlib library
 * published at www.zlib.net
 */


#define _FILE_OFFSET_BITS 64
#define _USE_KNETFILE


#include "bam.h"
#include "bgzf.h"

#include "blt_util/blt_exception.hh"
#include "boost/regex.hpp"

//#include <cassert>
//#include <clibio>
//#include <cctype>
#include <cstdio>
#include <cstdlib>
#include <iostream>

#include <string>


typedef int64_t __int64;

const unsigned MAX_CHROMOSOMES(50);

typedef struct {         // i = row
	int32_t genomic_pos; // 10001, 10001, 10001
	uint8_t nth;         // 0, 1, 2
} map_t;

typedef struct {
	char name[10]; //for marshalling fixed size helps
	int tid; //target id
	int start_index;
	int end_index;
	int n_ref;
	int genomic_pos_start;
	int genomic_pos_end;
	int n_alignments;
	int last_coverage_idx; //0 based
	bool has_data;
	int fd;
	float base_density; //coverage
	int approx_size_bytes;
	char *idl_path;
	bool dynamic;
	int err;
}chromosome_t;

typedef struct {
	BGZF *fp;
	const char *path;				//keep record of bam file path
	long n_alignments;		//number of alignments we are allowing
	map_t *maps;			//array of map_t
	bam_header_t *header;	//
	int fd;
	bam_index_t *loaded_index;
	uint32_t max_coverage_index;
	int real_max_linear_index;
	int *idx_value; 
	FILE *idxWrite; 
	FILE *idlWrite;
	FILE *idlRead;
	FILE *idxRead;
	int current_index_bound;
	int last_idx[25];
	bool dynamic;
	int avg_compressed_record_length;
	int record_visited;
	chromosome_t chromosomes[MAX_CHROMOSOMES];
	int active_tid;
} maps_t;


namespace GSBAMINDEX {

namespace BAM_SUBAPI {


#include "bam_index_subset.h"


static
bool bamw_has_data(const bam_index_t *idxt, int tid)
{
        const __bam_index_t* idx(reinterpret_cast<const __bam_index_t*>(idxt));
        khash_t(i) *index;
        index = idx->index[tid];
        if(index->n_buckets == 0 && index->size == 0)
                return false;
        return true;
}


static
bam1_t *bamw_fetch_one(bamFile fp, const bam_index_t *idxt, int tid, int beg, int end,  void * /*data*/, int genomic_pos, int nth, bool first)
{
        const __bam_index_t* idx(reinterpret_cast<const __bam_index_t*>(idxt));
	uint16_t *bins;
	int i, n_bins, n_off;
	pair64_t *off;
	khint_t k;
	khash_t(i) *index;
	uint64_t min_off;

	bins = (uint16_t*)calloc(BAM_MAX_BIN, 2);
	n_bins = reg2bins(beg, end, bins);
	index = idx->index[tid];
	min_off = (beg>>BAM_LIDX_SHIFT >= idx->index2[tid].n)? 0 : idx->index2[tid].offset[beg>>BAM_LIDX_SHIFT];
	for (i = n_off = 0; i < n_bins; ++i) {
		if ((k = kh_get(i, index, bins[i])) != kh_end(index))
			n_off += kh_value(index, k).n;
	}
	if (n_off == 0) {
		if(bins!= NULL) free(bins); return 0;
	}
	off = (pair64_t*)calloc(n_off, 16);
	for (i = n_off = 0; i < n_bins; ++i) {
		if ((k = kh_get(i, index, bins[i])) != kh_end(index)) {
			bam_binlist_t *p = &kh_value(index, k);
			for (unsigned j = 0; j < p->n; ++j)
				if (p->list[j].v > min_off) off[n_off++] = p->list[j];
		}
	}
	free(bins);
	{
		int positions = 0;
		bam1_t *b;
		int l, ret, n_seeks;
		uint64_t curr_off;
		b = (bam1_t*)calloc(1, sizeof(bam1_t));
		b->core.pos = -1;
		ks_introsort(off, n_off, off);
		// resolve completely contained adjacent blocks
		for (i = 1, l = 0; i < n_off; ++i)
			if (off[l].v < off[i].v)
				off[++l] = off[i];
		n_off = l + 1;
		// resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
		for (i = 1; i < n_off; ++i)
			if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
		{ // merge adjacent blocks
#if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16)
			for (i = 1, l = 0; i < n_off; ++i) {
#ifdef BAM_TRUE_OFFSET
				if (off[l].v + BAM_MIN_CHUNK_GAP > off[i].u) off[l].v = off[i].v;
#else
				if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
#endif
				else off[++l] = off[i];
			}
			n_off = l + 1;
#endif
		}
		// retrive alignments
		n_seeks = 0; i = -1; curr_off = 0;
		for (;;) {
			if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk
				if (i == n_off - 1) break; // no more chunks
				if (i >= 0) assert(curr_off == off[i].v); // otherwise bug
				if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek
					bam_seek(fp, off[i+1].u, SEEK_SET);
					curr_off = bam_tell(fp);
					++n_seeks;
				}
				++i;
			}
			if ((ret = bam_read1(fp, b)) > 0) {
				curr_off = bam_tell(fp);
				if (b->core.tid != tid || b->core.pos > end) break; // no need to proceed
				else if (is_overlap(beg, end, b)) 
				{	
					if(first && b->core.pos >= genomic_pos)
						return b;
						//return bam_format1((bam_header_t*)data, b);
					if(genomic_pos==b->core.pos && nth==positions)
						return b;
						//return bam_format1((bam_header_t*)data, b);
					positions++;
				}
			} else break; // end of file
		}
//		fprintf(stderr, "[bam_fetch] # seek calls: %d\n", n_seeks);
		if(first)
			return b; //last alignment
		bam_destroy1(b);
	}
	free(off);
	return 0;
}

} // namespace BAM_SUBAPI



const size_t map_t_size(sizeof(map_t));



/*    
   MA : decl should go in a header
   Any additional function we want to add as wrapper around original bam API, reside here
*/
const int MAX_BAM_FILES = 300;

//chromosome_t chromosomes[MAX_CHROMOSOMES]; 
//int i = 0;
char bases[17] = "=ACMGRSVTWYHKDBN";

//START - TEST API ONLY
char *str_seq;  
//realloc secret:
//Now suppose you've allocated a certain number of bytes for an array but later find that you want to add 
//values to it. You could copy everything into a larger array, which is inefficient, 
//or you can allocate more bytes using realloc, without losing your data.
int32_t size_starts = 1000000; //.1 million //TEST API ONLY
int32_t *starts = (int*)calloc(size_starts, sizeof(int32_t)); //calloc intializes to 0, malloc doesnt //TEST API ONLY
//END - TEST API ONLY

int32_t last_pos = 0;
const bam1_t *current_record;
uint8_t nth = 0;
static int32_t n = 0;
map_t user_map;

int indexBoundIncrement = 1000;    
maps_t arr_maps[MAX_BAM_FILES]; 

int fd_index; 
int current_fd_index = 0;


static int index_of(int fd)
{
	for(int i = 0; i < MAX_BAM_FILES; i++)
		if(arr_maps[i].fd == fd)
			return i;
	return -1;
}

static int index_of_chromosome(const char* name, const int fd)
{
	//
	maps_t *this_maps = &(arr_maps[index_of(fd)]);
	for(unsigned i = 0; i < MAX_CHROMOSOMES; i++)
		if(strcmp(name, this_maps->chromosomes[i].name)==0)
			return i;
	return -1;
}

//rewrite binary search
static int map_binary_search(int genomicPos, int lowIdx, int highIdx, int index_of_maps)
{	
	if(highIdx < lowIdx)
		return -1;
	int mid = lowIdx + (highIdx - lowIdx)/2;
	maps_t *this_maps = &(arr_maps[index_of_maps]);
	map_t *mid_user_map = (map_t*)calloc(1, sizeof(map_t));
	__int64 seekPos = (__int64)mid * sizeof(map_t);
	fseeko(this_maps->idlRead, seekPos, 0);
	const size_t n_read = fread(mid_user_map, map_t_size, 1, this_maps->idlRead);
        if(1 != n_read) {
            throw blt_exception("ERROR: Failed idl file read.\n");
        }
	if(mid_user_map->genomic_pos > genomicPos && mid > lowIdx)	
		return map_binary_search(genomicPos, lowIdx, mid-1, index_of_maps);
	else if(mid_user_map->genomic_pos < genomicPos && mid < highIdx)
		return map_binary_search(genomicPos, mid+1, highIdx, index_of_maps);
	else
		return mid;	
}

map_t common_map;
int last_tid;
int progressAt = 1000;
static int build_linear_index(const bam1_t *b, void * /*data*/)
{

    size_t n_write;

	maps_t *this_maps = &(arr_maps[current_fd_index]);
	if(n > 0 && common_map.genomic_pos == b->core.pos) // common_map.genomic_pos == last_pos
		nth++;
	else
		nth = 0;
	if(last_tid!= b->core.tid)	
	{
		if(last_tid > -1 && last_tid < this_maps->header->n_targets)
		{
			this_maps->chromosomes[last_tid].end_index = n-1;	
			this_maps->chromosomes[last_tid].genomic_pos_end = common_map.genomic_pos;
			this_maps->idx_value[this_maps->real_max_linear_index] = n;
			n_write = fwrite(&n, sizeof(int), 1, this_maps->idxWrite);
                        if(1 != n_write) {
                            throw blt_exception("ERROR: Failed idx file write.\n");
                        }
			this_maps->chromosomes[last_tid].last_coverage_idx = this_maps->real_max_linear_index;
			//this_maps->real_max_linear_index++;
		}
		 
		if(b->core.tid < this_maps->header->n_targets)
		{
			this_maps->chromosomes[b->core.tid].genomic_pos_start = b->core.pos;
			this_maps->chromosomes[b->core.tid].start_index = n;
			this_maps->current_index_bound = 0; // reset to zero for new chromosome
			this_maps->real_max_linear_index = 0;
		}
	}
	last_tid = b->core.tid;
	common_map.genomic_pos = b->core.pos;
	common_map.nth = nth;
	n_write = fwrite(&common_map, map_t_size, 1, this_maps->idlWrite);
        if(1 != n_write) {
            throw blt_exception("ERROR: Failed idl file write.\n");
        }
	while(b->core.pos > this_maps->current_index_bound)
	{
		//maps_t this_maps = arr_maps[current_fd_index]; // is it copying?
		this_maps->idx_value[this_maps->real_max_linear_index] = n;
		n_write = fwrite(&n, sizeof(int), 1, this_maps->idxWrite);
                if(1 != n_write) {
                    throw blt_exception("ERROR: Failed idx file write.\n");
                }
        this_maps->current_index_bound += indexBoundIncrement;
		if(b->core.tid == this_maps->header->n_targets -1 )
			this_maps->chromosomes[last_tid].last_coverage_idx = this_maps->real_max_linear_index;
		this_maps->real_max_linear_index++;		
	}
#if 0
	//progress per fd
	if(n==progressAt)
	{
		BGZF *f = this_maps->fp;
		__int64 cur_pos = ftello(f->file);	
		user_function_report_current_bamfile_pos(cur_pos, b->core.tid);
		progressAt += 1000;
	}
#endif
	n++;		
	return 0; 
}

#if 0
static
int bamw_close_active(int fd)
{
	maps_t *this_maps = &(arr_maps[current_fd_index]);
	if(this_maps->idlRead != NULL)
		fclose(this_maps->idlRead);
	if(this_maps->idxRead != NULL)
		fclose(this_maps->idxRead);
	return 0;
}
#endif


/*
with and without linear index
*/
static
int bamw_open_build_header(const char *bam_file_path, bool load_bai) //load_bai is almost always true
{
	const char mode[] = "rb"; //windows default isnt b, explicitly mention it as b
	maps_t *this_maps = &(arr_maps[fd_index]);
	if(load_bai)
		this_maps->loaded_index = bam_index_load(bam_file_path);
	this_maps->fp = bam_open(bam_file_path, mode);
	this_maps->header = bam_header_read(this_maps->fp);

        assert(this_maps->header->n_targets<=MAX_CHROMOSOMES);

#ifdef _USE_KNETFILE
	this_maps->fd = knet_fileno(this_maps->fp->x.fpr);
#else
	this_maps->fd = fileno(this_maps->fp->file);
#endif
	this_maps->path = bam_file_path;
	fd_index++;
	if(this_maps->header != NULL)
		for(int i=0; i< this_maps->header->n_targets; i++)
		{			
			//if(!chromosomes[i].has_data) //only consider those that has never been set to true, since we only a global set of chrs for all bam files
			{
				char *name= this_maps->header->target_name[i];				
				int ref = 0; int beg=0; int end = 0;					
				bam_parse_region(this_maps->header, name, &ref, &beg, &end); //ref == tid?
				strcpy(this_maps->chromosomes[i].name, name);
				this_maps->chromosomes[ref].n_ref = this_maps->header->target_len[i];
				this_maps->chromosomes[ref].tid = ref;
				this_maps->chromosomes[ref].has_data = false;
				if(load_bai) 							
                                    this_maps->chromosomes[ref].has_data = BAM_SUBAPI::bamw_has_data(this_maps->loaded_index, this_maps->chromosomes[ref].tid);
			}
			
		}

        return this_maps->fd;
}

/*
with and without linear index
*/
static
int bamw_init(const int fd, const char *index_file_path, bool create_ilmn_index, bool dynamic)
{
	last_tid = -1;
	current_fd_index = index_of(fd);
	maps_t *this_maps = &(arr_maps[current_fd_index]);
	
	int len = strlen(index_file_path);  
	for(int i=0; i<this_maps->header->n_targets ; i++)
		this_maps->max_coverage_index += (this_maps->header->target_len[i]) / indexBoundIncrement + 1;
	this_maps->idx_value = (int*)calloc(this_maps->max_coverage_index, sizeof(int));
	this_maps->dynamic = dynamic;
	if(!dynamic)
	{		
		char *idl_file_path = (char*)calloc(len, 1);
		strcpy(idl_file_path, index_file_path);
		idl_file_path[strlen(idl_file_path)-1] = 'l';

		if(create_ilmn_index)
		{
			progressAt = 1000;
			this_maps->idxWrite = fopen(/*this_maps->path*/index_file_path,"wb"); //for whatever reason ->path is lost between fucntion calls
			this_maps->idlWrite = fopen(idl_file_path, "wb");
			int ref = 0; int beg=0; int end = 0;
			bam_fetch_f result_func = &build_linear_index;
			for(int j=0; j<this_maps->header->n_targets; j++)
			{				
				if(!this_maps->chromosomes[j].has_data)
					continue;
				bam_parse_region(this_maps->header, this_maps->chromosomes[j].name, &ref, &beg, &end); //ref == tid?
				this_maps->chromosomes[j].tid = ref;
				bam_fetch(this_maps->fp, this_maps->loaded_index, ref, beg, end, this_maps->header, result_func); //ecoli == 0?
				this_maps->chromosomes[j].dynamic = false;
			}	
			fclose(this_maps->idxWrite);
			fclose(this_maps->idlWrite);
			
			this_maps->chromosomes[last_tid].end_index = n;	
			this_maps->chromosomes[last_tid].genomic_pos_end = common_map.genomic_pos;

			this_maps->idlRead = fopen(idl_file_path, "rb");
			this_maps->idxRead = fopen(index_file_path, "rb");

			map_t *start_user_map = (map_t*)calloc(1, sizeof(map_t)); 
			__int64 seekPos = (__int64)0 * sizeof(map_t);
			fseeko(this_maps->idlRead, seekPos, 0);
			const size_t n_read = fread(start_user_map, map_t_size, 1, this_maps->idlRead);
                        if(1 != n_read) {
                            throw blt_exception("ERROR: Failed idl file read.\n");
                        }

			this_maps->chromosomes[0].genomic_pos_start = start_user_map->genomic_pos; //we never updated start genomic pos for tid 0 in the call back 
			this_maps->n_alignments = n;
			//chromosomes[this_maps->header->n_targets-1].last_coverage_idx = this_maps->real_max_linear_index;
			n = 0;

                        free(start_user_map);
		}
		else
		{
			this_maps->idlRead = fopen(idl_file_path, "rb");
			this_maps->idxRead = fopen(index_file_path, "rb");
		}
	}
	else //dynamic
	{
		//int ref = 0; int beg=0; int end = 0;
		for(int j=0; j<this_maps->header->n_targets; j++)
		{
			if(!this_maps->chromosomes[j].has_data)
				continue;
			int ref = 0; int beg=0; int end = 0;
			bam_parse_region(this_maps->header, this_maps->chromosomes[j].name, &ref, &beg, &end); //ref == tid?
			this_maps->chromosomes[j].tid = ref;
			bam1_t *first = BAM_SUBAPI::bamw_fetch_one(this_maps->fp, this_maps->loaded_index, ref, beg, end, this_maps->header, user_map.genomic_pos, user_map.nth, true);
			//need better measure, typically 1 compressed byte for each read base. 
			//int size1 = first->data_len + sizeof(bam1_core_t) + 3*sizeof(int);
			int size1 = first->core.l_qseq;
			this_maps->chromosomes[j].genomic_pos_end = this_maps->header->target_len[ref];
			int last_pos = this_maps->chromosomes[j].genomic_pos_end;
			if(first != NULL && first->core.pos > -1)
				this_maps->chromosomes[j].genomic_pos_start = first->core.pos;
			else
			{				
				this_maps->chromosomes[j].err = 1;
				continue; 
			}
			
			beg=this_maps->chromosomes[j].genomic_pos_end; 
			end = this_maps->chromosomes[j].genomic_pos_end;
			
			bam1_t *last =  BAM_SUBAPI::bamw_fetch_one(this_maps->fp, this_maps->loaded_index, ref, beg, end, this_maps->header, last_pos, 0, true);			
			this_maps->chromosomes[j].dynamic = true;
			int max = -1;
			if(last != NULL && last->core.pos > -1)
				max = last->core.pos;
			else
			{
				//try binary search like search!
				int start_search = first->core.pos;
				int end_search = last_pos;
				int mid_search = (start_search + end_search)/2;
				
				while(mid_search > start_search)
				{
                                    last =  BAM_SUBAPI::bamw_fetch_one(this_maps->fp, this_maps->loaded_index, ref, mid_search, end_search, this_maps->header, mid_search, 0, true);
					if(last == NULL || last->core.pos < mid_search)
					{
						end_search = mid_search;
						mid_search = (start_search + end_search)/2;
					}
					else
					{
						if(last->core.pos > max)
							max = last->core.pos;
						mid_search = last->core.pos;
						start_search = mid_search;
						mid_search = (start_search + end_search)/2;
					}
				}				
			}

			if(last==NULL || last->core.pos < 0)
                            last = BAM_SUBAPI::bamw_fetch_one(this_maps->fp, this_maps->loaded_index, ref, max, max+1000, this_maps->header, max, 0, true);
			this_maps->chromosomes[j].genomic_pos_end = last->core.pos;
			int size2 = last->core.l_qseq;
			this_maps->avg_compressed_record_length = (size1 + size2)/2;

			if(last != NULL)
				this_maps->chromosomes[j].genomic_pos_end = last->core.pos;
			free(last);
			last = NULL;			
			
			BGZF *f = this_maps->fp;
                        //__int64 cur_pos;
#ifdef _USE_KNETFILE                
                        knet_tell(f->x.fpr);
#else
			ftello(f->file);
#endif
		}	
	}
	//	
	//n = 0;
	int idx_fd = -1;
	if(this_maps->idlRead != NULL)
		idx_fd = fileno(this_maps->idxRead);
	return idx_fd;
}

/*
 client should maintain fd, should expect array of int with two fd,  one for idx, one for idl
*/
static
void bamw_create_linear_index(int bam_fd, const char *chr, const char *idx_file_path)
{
	last_tid = -1;
	current_fd_index = index_of(bam_fd);
	maps_t *this_maps = &(arr_maps[current_fd_index]);
	const int index = index_of_chromosome(chr, bam_fd);
	
	chromosome_t *chromosome = &(this_maps->chromosomes[index]);
        {
            const int len(strlen(idx_file_path));
            chromosome->idl_path = (char*)calloc((len+1), 1);
            strcpy(chromosome->idl_path, idx_file_path);
            chromosome->idl_path[len-1] = 'l';
        }
        const char* idl_file_path(chromosome->idl_path);
	progressAt = 1000;
	this_maps->idxWrite = fopen(idx_file_path,"wb"); 
	this_maps->idlWrite = fopen(idl_file_path, "wb");
	this_maps->active_tid = chromosome->tid;
	//this_maps->idxRead = fopen(index_file_path, "rb");
	int ref = 0; int beg=0; int end = 0;
	bam_fetch_f result_func = &build_linear_index;

	bam_parse_region(this_maps->header, chr, &ref, &beg, &end);
	chromosome->tid = ref;
	
	bam_fetch(this_maps->fp, this_maps->loaded_index, ref, beg, end, this_maps->header, result_func); 
	fclose(this_maps->idxWrite);
	fclose(this_maps->idlWrite);
	
	chromosome->end_index = n;	
	chromosome->genomic_pos_end = common_map.genomic_pos;
	chromosome->dynamic = false; 

	this_maps->idlRead = fopen(idl_file_path, "rb");
	map_t *start_user_map = (map_t*)calloc(1, sizeof(map_t)); 
	__int64 seekPos = (__int64)0 * sizeof(map_t);
	fseeko(this_maps->idlRead, seekPos, 0);
	const size_t n_read = fread(start_user_map, map_t_size, 1, this_maps->idlRead);
        if(1 != n_read) {
            throw blt_exception("ERROR: Failed idl file read.\n");
        }

	chromosome->genomic_pos_start = start_user_map->genomic_pos; //we never updated start genomic pos for tid 0 in the call back 
	this_maps->n_alignments += n;

        free(start_user_map);
	fclose(this_maps->idlRead);
	n = 0;
}

/*
with and without linear index
*/
static
chromosome_t *bamw_get_target_info(const char *chromosome, int fd)
{
	maps_t *this_maps = &(arr_maps[index_of(fd)]);
	const int index = index_of_chromosome(chromosome, fd);
	return &this_maps->chromosomes[index];
}

#if 0
/*
with and without linear index
*/
static
const char *bamw_get_header_text(int fd)
{
	if(arr_maps[index_of(fd)].header!=NULL)
		return arr_maps[index_of(fd)].header->text;
	return NULL;
}
#endif

/*
with and without linear index
*/
static
bam_header_t *bamw_get_header(int fd)
{
	return (arr_maps[index_of(fd)].header);
}

#if 0
/*
with and without linear index
*/
static
bam_header_t *bamw_get_header_close(const char *bam_file_path) //and close
{
	const char mode[] = "rb"; //windows default isnt b, explicitly mention it as b
	BGZF *fp = bam_open(bam_file_path, mode);
	bam_header_t *header = bam_header_read(fp);
	bam_index_t *bai_index = bam_index_load(bam_file_path);
	
	if(header != NULL)
		for(int i=0; i< header->n_targets; i++)
		{		
			char *name= header->target_name[i];
			int ref = 0; int beg=0; int end = 0;
			bam_parse_region(header, name, &ref, &beg, &end); //ref == tid?
			if(!BAM_SUBAPI::bamw_has_data(bai_index, ref))			
				header->target_name[i]= NULL;			
		}
	bam_index_destroy(bai_index);
	//validate targets
	bam_close(fp);
	return header;
}
#endif

/*
with and without linear index
*/
static
char **bamw_get_targets(int fd)
{
	if(arr_maps[index_of(fd)].header != NULL)
		return arr_maps[index_of(fd)].header->target_name;
	return NULL;
}

} // namespace GSBAMINDEX


static
const char* change_file_ext(const char* file_path, const char* ext)
{
	size_t file_path_len = strlen(file_path);
	char* new_file_path = (char*) calloc(file_path_len+1, 1);
	strncpy(new_file_path, file_path, file_path_len-3); 
	new_file_path[file_path_len-3]='\0';
	strcat(new_file_path, ext);
	return new_file_path;
}

static
void normalize_chrom_name(std::string& chrom){

    static const boost::regex reLead("^c(hr)?");//, boost::regex_constants::icase);
    static const boost::regex reTrail(".fa$");
    static const boost::regex reMito("^mi?t?o?c?h?o?n?d?r?i?((on)|(al)|(a))?");

    const unsigned cs(chrom.size());
    for(unsigned i(0);i<cs;++i) chrom[i] = tolower(chrom[i]);
  
    chrom=regex_replace(chrom, reLead,"");
    chrom=regex_replace(chrom, reTrail,"");
    chrom=regex_replace(chrom, reMito,"m");
}

static
void xml_serialize_start(const char* bam_file_path,FILE*& xml_file)
{
        const char* xml_file_path = change_file_ext(bam_file_path, "xml");
        xml_file = fopen(xml_file_path, "w");
        fprintf(xml_file, "%s", "<?xml version=\"1.0\"?>\n");
        fprintf(xml_file, "%s", "<ArrayOfTarget_t xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xmlns:xsd=\"http://www.w3.org/2001/XMLSchema\">\n");
}

static
void xml_serialize_end(FILE* xml_file)
{
	fprintf(xml_file, "%s", "</ArrayOfTarget_t>\n");
	fclose(xml_file);
}

static
void xml_serialize_target(const chromosome_t* t,FILE* xml_file)
{
	fprintf(xml_file, "%s", "  <target_t>\n");
	fprintf(xml_file, "    <name>%s</name>\n",t->name);
	fprintf(xml_file, "    <tid>%u</tid>\n",t->tid);
	fprintf(xml_file, "    <start_index>%u</start_index>\n",t->start_index);
	fprintf(xml_file, "    <end_index>%u</end_index>\n",t->end_index);
	fprintf(xml_file, "    <n_ref>%u</n_ref>\n",t->n_ref);
	fprintf(xml_file, "    <genomic_pos_start>%u</genomic_pos_start>\n",t->genomic_pos_start);
	fprintf(xml_file, "    <genomic_pos_end>%u</genomic_pos_end>\n",t->genomic_pos_end);
	fprintf(xml_file, "    <n_alignments>%u</n_alignments>\n",t->n_alignments);
	fprintf(xml_file, "    <last_coverage_idx>%u</last_coverage_idx>\n",t->last_coverage_idx);
	fprintf(xml_file, "    <has_data>%s</has_data>\n",t->has_data?"true":"false");
	fprintf(xml_file, "    <fd>%u</fd>\n",t->fd);
	fprintf(xml_file, "    <baseDensity>%g</baseDensity>\n",t->base_density);
	fprintf(xml_file, "    <approx_size_bytes>%u</approx_size_bytes>\n",t->approx_size_bytes);
	fprintf(xml_file, "    <dynamic>%s</dynamic>\n",t->dynamic?"true":"false");
	fprintf(xml_file, "    <err>%i</err>\n",t->err);
	fprintf(xml_file, "%s", "  </target_t>\n");
}


static
void
try_main(int argc, char* argv[])
{
        const char* progname(argv[0]);

        if(argc!=2) {
            fprintf(stderr,"\nusage: %s bamfile\n\n",progname);
            exit(EXIT_FAILURE);
        }

        const char* bam_file_path(argv[1]);

	const char* index_file_path = change_file_ext(bam_file_path, "idx");

        std::string bam_file_dir;
        {
            const char* x(strrchr(bam_file_path,'/'));
            if(NULL!=x) {
                bam_file_dir = std::string(bam_file_path,x+1);
            }
        }

	int bam_fd = GSBAMINDEX::bamw_open_build_header(bam_file_path, true);
	
        GSBAMINDEX::bamw_init(bam_fd, index_file_path, false, true);

	bam_header_t* bam_header = GSBAMINDEX::bamw_get_header(bam_fd);

        FILE* xml_file;
	xml_serialize_start(bam_file_path,xml_file);

	// getting a list of chromosomes
	char** chromosomes = GSBAMINDEX::bamw_get_targets(bam_fd);

	// creating an index file for each chromosome
	for(int i=0; i<bam_header->n_targets; i++)
	{
            assert(NULL != chromosomes);
            assert(NULL != chromosomes[i]);

		std::string chrom(chromosomes[i]);
		chromosome_t* target_info = GSBAMINDEX::bamw_get_target_info(chrom.c_str(), bam_fd);

                normalize_chrom_name(chrom);
                const std::string idx_file_path(bam_file_dir+chrom+".idx");
                if(target_info->has_data) {
                    GSBAMINDEX::bamw_create_linear_index(bam_fd, target_info->name, idx_file_path.c_str());
                    xml_serialize_target(target_info,xml_file);
                }

	}

	xml_serialize_end(xml_file);
}



static
void
dump_cl(int argc,
        char* argv[],
        std::ostream& os) {

    os << "cmdline:";
    for(int i(0);i<argc;++i){
        os << ' ' << argv[i];
    }
    os << std::endl;
}



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

    const char* progname(argv[0]);

    std::ostream& log_os(std::cerr);

    // last chance to catch exceptions...
    //
    try{
        try_main(argc,argv);

    } catch (const blt_exception& e) {
        log_os << "FATAL_ERROR: " << progname << " EXCEPTION: " << e.what() << "\n"
               << "...caught in main()\n";
        dump_cl(argc,argv,log_os);
        exit(EXIT_FAILURE);

    } catch(const std::exception& e) {
        log_os << "FATAL_ERROR: EXCEPTION: " << e.what() << "\n"
               << "...caught in main()\n";
        dump_cl(argc,argv,log_os);
        exit(EXIT_FAILURE);

    } catch(...) {
        log_os << "FATAL_ERROR: UNKNOWN EXCEPTION\n"
               << "...caught in main()\n";
        dump_cl(argc,argv,log_os);
        exit(EXIT_FAILURE);
    }
    return EXIT_SUCCESS;
}