#ifdef PERL_CAPI
#define WIN32IO_IS_STDIO
#endif
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#ifdef FCGI
 #include <fcgi_stdio.h>
#else
 #ifdef USE_SFIO
  #include <config.h>
 #else
  #include <stdio.h>
 #endif
 #include <perlio.h>
#endif

#ifndef Newx
#  define Newx(v,n,t) New(0,v,n,t)
#endif

#ifndef Newxz
#  define Newxz(v,n,t) Newz(0,v,n,t)
#endif

#include <unistd.h>
#include <math.h>
#include "bam.h"
#include "khash.h"
#include "faidx.h"

/* stolen from bam_aux.c */
#define MAX_REGION 1<<29

/* stolen from bam.c */
extern void bam_view1(const bam_header_t *header, const bam1_t *b);
/* stolen from bam_sort.c */
extern void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t max_mem);

typedef tamFile         Bio__DB__Tam;
typedef faidx_t*        Bio__DB__Sam__Fai;
typedef bamFile         Bio__DB__Bam;
typedef bam_header_t*   Bio__DB__Bam__Header;
typedef bam1_t*         Bio__DB__Bam__Alignment;
typedef bam_index_t*    Bio__DB__Bam__Index;
typedef bam_pileup1_t*  Bio__DB__Bam__Pileup;
typedef struct {
  SV* callback;
  SV* data;
} fetch_callback_data;
typedef fetch_callback_data *fetch_callback_dataptr;
typedef struct {
  int    start;
  int    end;
  double width;
  int    reads;
  int*   bin;
} coverage_graph;
typedef coverage_graph *coverage_graph_ptr;

static int MaxPileupCnt=8000;

void XS_pack_charPtrPtr( SV * arg, char ** array, int count) {
  int i;
  AV * avref;
  avref = (AV*)sv_2mortal((SV*)newAV());
  for (i=0; i<count; i++) {
    av_push(avref, newSVpv(array[i], strlen(array[i])));
  }
  SvSetSV( arg, newRV((SV*)avref));
}

int bam_fetch_fun (const bam1_t *b, void *data) {
  dSP;
  int count;

  fetch_callback_dataptr fcp;
  SV* callback;
  SV* callbackdata;
  SV* alignment_obj;
  bam1_t *b2;

  fcp          = (fetch_callback_dataptr) data;
  callback     = fcp->callback;
  callbackdata = fcp->data;

  /* turn the bam1_t into an appropriate object */
  /* need to dup it here so that the C layer doesn't reuse the address under Perl */
  b2 = bam_dup1(b);

  alignment_obj = sv_setref_pv(newSV(sizeof(bam1_t)),"Bio::DB::Bam::Alignment",(void*) b2);

  /* set up subroutine stack for the call */
  ENTER;
  SAVETMPS;

  PUSHMARK(SP);
  XPUSHs(sv_2mortal(alignment_obj));
  XPUSHs(callbackdata);
  PUTBACK;

  /* execute the call */
  count = call_sv(callback,G_SCALAR|G_DISCARD);

  FREETMPS;
  LEAVE;

  return 1;
}

int invoke_pileup_callback_fun(uint32_t tid, 
			       uint32_t pos, 
			       int n, 
			       const bam_pileup1_t *pl,
			       void *data) {
  dSP;
  int count,i;
  fetch_callback_dataptr fcp;
  SV*  callback;
  SV*  callbackdata;
  SV*  pileup_obj;
  SV* p;
  SV** pileups;
  AV*  pileup;

  fcp          = (fetch_callback_dataptr) data;
  callback     = fcp->callback;
  callbackdata = fcp->data;

  /* turn the bam_pileup1_t into the appropriate object */
  /* this causes a compiler warning -- ignore it */
  pileup = newAV();
  av_extend(pileup,n);
  for (i=0;i<n;i++) {
    p = newSV(sizeof(bam_pileup1_t));	
    sv_setref_pv(p,"Bio::DB::Bam::Pileup",(void*) &pl[i]);
    av_push(pileup,p);
  } 
  
  /* set up subroutine stack for the call */
  ENTER;
  SAVETMPS;

  PUSHMARK(SP);
  XPUSHs(sv_2mortal(newSViv(tid)));
  XPUSHs(sv_2mortal(newSViv(pos)));
  XPUSHs(sv_2mortal(newRV_noinc((SV*)pileup)));
  XPUSHs(callbackdata);
  PUTBACK;

  /* execute the call */
  count = call_sv(callback,G_SCALAR|G_DISCARD);

  FREETMPS;
  LEAVE;
}

int add_pileup_line (const bam1_t *b, void *data) {
  bam_plbuf_t *pileup = (bam_plbuf_t*) data;
  bam_plbuf_push(b,pileup);
  return 0;
}

int add_lpileup_line (const bam1_t *b, void *data) {
  bam_lplbuf_t *pileup = (bam_lplbuf_t*) data;
  bam_lplbuf_push(b,pileup);
  return 0;
}

int coverage_from_pileup_fun (uint32_t tid, 
			      uint32_t pos, 
			      int n, 
			      const bam_pileup1_t *pl, 
			      void *data) {
  coverage_graph_ptr  cgp;
  int                 bin;
  int                 i;
  int                 valid;

  cgp = (coverage_graph_ptr) data;
  cgp->reads += n;

  valid = 0;
  for (i=0;i<n;i++) {
    if (!pl[i].is_del && !pl[i].is_refskip)
        valid++;
  }

  if (pos >= cgp->start && pos <= cgp->end) {
    bin = (pos-cgp->start)/cgp->width;
    cgp->bin[bin] += valid;
  }

  return 0;
}

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Tam PREFIX=tam_

Bio::DB::Tam
tam_open(packname="Bio::DB::Tam", filename)
   char * packname
   char * filename
 PROTOTYPE: $$
 CODE:
    RETVAL = sam_open(filename);
 OUTPUT:
    RETVAL

void
tam_DESTROY(tam)
  Bio::DB::Tam tam
  PROTOTYPE: $
  CODE:
     sam_close(tam);

Bio::DB::Bam::Header
tam_header_read2(packname="Bio::DB::Tam", filename)
    char * packname
    char * filename
    PROTOTYPE: $$
    CODE:
      RETVAL = sam_header_read2(filename);
    OUTPUT:
      RETVAL

Bio::DB::Bam::Header
tam_header_read(tam)
    Bio::DB::Tam            tam
    PROTOTYPE: $$
    CODE:
      RETVAL = sam_header_read(tam);
    OUTPUT:
      RETVAL

int
tam_read1(tam,header,alignment)
    Bio::DB::Tam            tam
    Bio::DB::Bam::Header    header
    Bio::DB::Bam::Alignment alignment
    CODE:
       RETVAL = sam_read1(tam,header,alignment);
    OUTPUT:
       RETVAL

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Sam::Fai PREFIX=fai_

Bio::DB::Sam::Fai
fai_load(packname="Bio::DB::Sam::Fai", filename)
  char * packname
  char * filename
 PROTOTYPE: $$
 CODE:
    RETVAL = fai_load(filename);
 OUTPUT:
    RETVAL

void
fai_destroy(fai)
  Bio::DB::Sam::Fai fai
  PROTOTYPE: $
  CODE:
    fai_destroy(fai);

SV*
fai_fetch(fai,reg)
  Bio::DB::Sam::Fai    fai
    const char *reg
  PROTOTYPE: $$$
  PREINIT:
    char     *seq;
    int       len;
  CODE:
    seq = fai_fetch(fai,reg,&len);
    if (seq == NULL)
       XSRETURN_EMPTY;
    RETVAL = newSVpv(seq,len);
    free((void*)seq);
  OUTPUT:
    RETVAL


MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Bam PREFIX=bam_

int
max_pileup_cnt(packname,...)
CODE:
	if (items > 1)
	   MaxPileupCnt = SvIV(ST(1));
	RETVAL = MaxPileupCnt;
OUTPUT:
        RETVAL

Bio::DB::Bam
bam_open(packname, filename, mode="r")
      char * packname
      char * filename
      char * mode
      PROTOTYPE: $$$
      CODE:
        RETVAL = bam_open(filename,mode);
      OUTPUT:
      RETVAL

void
bam_DESTROY(bam)
   Bio::DB::Bam bam
PROTOTYPE: $
CODE:
   bam_close(bam);

int
bam_index_build(packname, filename)
   char *      packname
   const char * filename
  CODE:
     RETVAL = bam_index_build(filename);
  OUTPUT:
     RETVAL

void
bam_sort_core(packname, is_by_qname=0, filename, prefix, max_mem=500000000)
   char * packname
   int    is_by_qname
   char * filename
   char * prefix
   int    max_mem
 PROTOTYPE: $$$$$
 CODE:
   bam_sort_core(is_by_qname,filename,prefix,max_mem);

Bio::DB::Bam::Index
bam_index_open(packname="Bio::DB::Bam", filename)
      char * packname
      char * filename
    PROTOTYPE: $$
    CODE:
    RETVAL = bam_index_load(filename);
    OUTPUT:
    RETVAL

Bio::DB::Bam::Header
bam_header(bam)
    Bio::DB::Bam bam
    PROTOTYPE: $
    PREINIT:
      bam_header_t *bh;
      int64_t       result;
    CODE:
      result = bgzf_seek(bam,0,0);
      bh = bam_header_read(bam);
      RETVAL = bh;
    OUTPUT:
      RETVAL

int
bam_header_write(bam,header)
    Bio::DB::Bam         bam
    Bio::DB::Bam::Header header
    PROTOTYPE: $$
    CODE:
      bgzf_seek(bam,0,0);
      RETVAL= bam_header_write(bam,header);
    OUTPUT:
      RETVAL

char*
bam_tell(bam)
    Bio::DB::Bam bam
PROTOTYPE: $
CODE:
    int64_t t = bam_tell(bam);
    char    string[128];
sprintf(string,"%llu",(long long unsigned int) t);
    RETVAL = string;
OUTPUT:
    RETVAL

void
bam_seek(bam,pos,dir)
    Bio::DB::Bam bam
    int pos
    int dir
PROTOTYPE: $$$
CODE:
    bam_seek(bam,pos,dir);

Bio::DB::Bam::Alignment
bam_read1(bam)
    Bio::DB::Bam bam
  PROTOTYPE: $
  PREINIT:
    bam1_t *b;
  CODE:
    b = bam_init1();
    if (bam_read1(bam,b) >= 0) {
      RETVAL = b;
    }
    else
       XSRETURN_EMPTY;
  OUTPUT:
    RETVAL      

int
bam_write1(bam,align)
    Bio::DB::Bam            bam
    Bio::DB::Bam::Alignment align
   PROTOTYPE: $$
   CODE:
      RETVAL = bam_write1(bam,align);
   OUTPUT:
      RETVAL

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Bam::Alignment PREFIX=bama_

Bio::DB::Bam::Alignment
bama_new(package="Bio::DB::Bam::Alignment")
   char * package
   PROTOTYPE: $
   CODE:
      RETVAL = bam_init1();
   OUTPUT:
      RETVAL

void
bama_DESTROY(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    bam_destroy1(b);

int
bama_tid(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.tid = SvIV(ST(1));
    RETVAL=b->core.tid;
OUTPUT:
    RETVAL

int
bama_pos(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.pos = SvIV(ST(1));
    RETVAL=b->core.pos;
OUTPUT:
    RETVAL

int
bama_calend(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
   RETVAL=bam_calend(&b->core,bam1_cigar(b));
OUTPUT:
   RETVAL    

int
bama_cigar2qlen(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
   RETVAL=bam_cigar2qlen(&b->core,bam1_cigar(b));
OUTPUT:
   RETVAL    

int
bama_qual(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.qual = SvIV(ST(1));
    RETVAL=b->core.qual;
OUTPUT:
    RETVAL

int
bama_flag(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.flag = SvIV(ST(1));
    RETVAL=b->core.flag;
OUTPUT:
    RETVAL

int
bama_n_cigar(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
  if (items > 1)
    b->core.n_cigar = SvIV(ST(1));
    RETVAL=b->core.n_cigar;
OUTPUT:
    RETVAL

int
bama_l_qseq(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.l_qseq = SvIV(ST(1));
    RETVAL=b->core.l_qseq;
OUTPUT:
    RETVAL

SV*
bama_qseq(b)
Bio::DB::Bam::Alignment b
PROTOTYPE: $
PREINIT:
    char* seq;
    int   i;
CODE:
    seq = Newxz(seq,b->core.l_qseq+1,char);
    for (i=0;i<b->core.l_qseq;i++) {
      seq[i]=bam_nt16_rev_table[bam1_seqi(bam1_seq(b),i)];
    }
    RETVAL = newSVpv(seq,b->core.l_qseq);
    Safefree(seq);
OUTPUT:
    RETVAL

SV*
bama__qscore(b)
Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL = newSVpv(bam1_qual(b),b->core.l_qseq);
OUTPUT:
    RETVAL

int
bama_mtid(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.mtid = SvIV(ST(1));
    RETVAL=b->core.mtid;
OUTPUT:
    RETVAL

int
bama_mpos(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.mpos = SvIV(ST(1));
    RETVAL=b->core.mpos;
OUTPUT:
    RETVAL

int
bama_isize(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.isize = SvIV(ST(1));
    RETVAL=b->core.isize;
OUTPUT:
    RETVAL

int
bama_l_aux(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->l_aux = SvIV(ST(1));
    RETVAL=b->l_aux;
OUTPUT:
    RETVAL

char*
bama_aux(b)
   Bio::DB::Bam::Alignment b
PREINIT:
   uint8_t *s;
   uint8_t type, key[2];
   char    str[8192];
CODE:
   s = bam1_aux(b);
   str[0] = '\0';

   int left  = sizeof(str) - strlen(str);
   while (left > 0 && (s < b->data + b->data_len)) {
        char* d   = str+strlen(str); 

	key[0] = s[0]; 
	key[1] = s[1];
 	left -= snprintf(d, left, "%c%c:", key[0], key[1]);

	d    += 3;
	s    += 2;
	type = *s++; 

	if (left <= 0) continue;

	if (type == 'A')      { left -= snprintf(d, left, "A:%c", *s);           s++; }
	else if (type == 'C') { left -= snprintf(d, left, "i:%u", *s);           s++; }
	else if (type == 'c') { left -= snprintf(d, left, "i:%d", *s);           s++; }
	else if (type == 'S') { left -= snprintf(d, left, "i:%u", *(uint16_t*)s);s += 2; }
	else if (type == 's') { left -= snprintf(d, left, "i:%d", *(int16_t*)s); s += 2; }
	else if (type == 'I') { left -= snprintf(d, left, "i:%u", *(uint32_t*)s);s += 4; }
	else if (type == 'i') { left -= snprintf(d, left, "i:%d", *(int32_t*)s); s += 4; }
	else if (type == 'f') { left -= snprintf(d, left, "f:%g", *(float*)s);   s += 4; }
	else if (type == 'd') { left -= snprintf(d, left, "d:%lg", *(double*)s); s += 8; }
	else if (type == 'Z' || type == 'H') { left -= snprintf(d, left, "%c:", type); 
	                                       strncat(d,s,left);
					       while (*s++) {}
					       left = sizeof(str) - strlen(str);
	                                     }
	if (left <= 0) continue;
	strncat(d,"\t",left);
	left--;
   }	  
   str[strlen(str)-1] = '\0';
   RETVAL = str;
OUTPUT:
   RETVAL

SV*
bama_aux_get(b,tag)
   Bio::DB::Bam::Alignment b
   char*               tag
PROTOTYPE: $$
PREINIT:
   int           type;
   uint8_t       *s;
CODE:
   s    = bam_aux_get_core(b,tag);
   if (s==0)
      XSRETURN_EMPTY;
   type = *s++;
   switch (type) {
   case 'c':
     RETVAL = newSViv((int32_t)*(int8_t*)s);
     break;
   case 'C':
     RETVAL = newSViv((int32_t)*(uint8_t*)s);
     break;
   case 's':
     RETVAL = newSViv((int32_t)*(int16_t*)s);
     break;
   case 'S':
     RETVAL = newSViv((int32_t)*(uint16_t*)s);
     break;
   case 'i':
     RETVAL = newSViv(*(int32_t*)s);
     break;
   case 'I':
     RETVAL = newSViv((int32_t)*(uint32_t*)s);
     break;
   case 'f':
     RETVAL = newSVnv(*(float*)s);
     break;
   case 'Z':
   case 'H':
     RETVAL = newSVpv((char*)s,0);
     break;
   case 'A':
     RETVAL = newSVpv((char*)s,1);
     break;
   default:
     XSRETURN_EMPTY;
   }
OUTPUT:
   RETVAL

void
bama_aux_keys(b)
Bio::DB::Bam::Alignment b
PROTOTYPE: $
PREINIT:
   uint8_t *s;
   uint8_t type;
PPCODE:
   {
     s = bam1_aux(b);  /* s is a khash macro */
     while (s < b->data + b->data_len) {
       XPUSHs(sv_2mortal(newSVpv(s,2)));
       s   += 2; 
       type = *s++;
       if      (type == 'A') { ++s; }
       else if (type == 'C') { ++s; }
       else if (type == 'c') { ++s; }
       else if (type == 'S') { s += 2; }
       else if (type == 's') { s += 2; }
       else if (type == 'I') { s += 4; }
       else if (type == 'i') { s += 4; }
       else if (type == 'f') { s += 4; }
       else if (type == 'Z' || type == 'H') { while (*s) ++(s); ++(s); }
     }
   }

SV*
bama_data(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
PREINIT:
    STRLEN  len;
CODE:
    if (items > 1) {
      b->data     = SvPV(ST(1),len);
      b->data_len = len;
    }
    RETVAL=newSVpv(b->data,b->data_len);
OUTPUT:
    RETVAL

int
bama_data_len(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->data_len = SvIV(ST(1));
    RETVAL=b->data_len;
OUTPUT:
    RETVAL

int
bama_m_data(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1) {
      b->m_data = SvIV(ST(1));
    }
    RETVAL=b->m_data;
OUTPUT:
    RETVAL

SV*
bama_qname(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL=newSVpv(bam1_qname(b),0);
OUTPUT:
    RETVAL

int
bama_paired(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL=(b->core.flag&BAM_FPAIRED) != 0;
OUTPUT:
  RETVAL

int
bama_proper_pair(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL=(b->core.flag&BAM_FPROPER_PAIR) != 0;
OUTPUT:
  RETVAL

int
bama_unmapped(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL=(b->core.flag&BAM_FUNMAP) != 0;
OUTPUT:
  RETVAL

int
bama_munmapped(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
    RETVAL=(b->core.flag&BAM_FMUNMAP) != 0;
OUTPUT:
  RETVAL

int
bama_reversed(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
  RETVAL=bam1_strand(b);
OUTPUT:
  RETVAL

int
bama_mreversed(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
CODE:
  RETVAL=bam1_mstrand(b);
OUTPUT:
  RETVAL

SV*
bama_cigar(b)
  Bio::DB::Bam::Alignment b
PROTOTYPE: $
PREINIT:
    int        i;
    uint32_t  *c;
    AV        *avref;
CODE:
    avref = (AV*) sv_2mortal((SV*)newAV());
    c     = bam1_cigar(b);
    for (i=0;i<b->core.n_cigar;i++)
      av_push(avref, newSViv(c[i]));
    RETVAL = (SV*) newRV((SV*)avref); 
OUTPUT:
  RETVAL

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Bam::Header PREFIX=bam_

Bio::DB::Bam::Header
bam_new(packname=Bio::DB::Bam::Header)
PROTOTYPE: $
CODE:
    RETVAL = bam_header_init();
OUTPUT:
    RETVAL

int
bam_n_targets(bamh)
  Bio::DB::Bam::Header bamh
  PROTOTYPE: $
  CODE:
    RETVAL = bamh->n_targets;
  OUTPUT:
    RETVAL

SV*
bam_target_name(bamh)
  Bio::DB::Bam::Header bamh
  PROTOTYPE: $
  PREINIT:
    int i;
    AV * avref;
  CODE:
    avref = (AV*) sv_2mortal((SV*)newAV());
    for (i=0;i<bamh->n_targets;i++)
      av_push(avref, newSVpv(bamh->target_name[i],0));
    RETVAL = (SV*) newRV((SV*)avref); 
  OUTPUT:
    RETVAL

SV*
bam_target_len(bamh)
    Bio::DB::Bam::Header bamh
  PROTOTYPE: $
  PREINIT:
    int i;
    AV * avref;
  CODE:
    avref = (AV*) sv_2mortal((SV*)newAV());
    for (i=0;i<bamh->n_targets;i++)
       av_push(avref, newSViv(bamh->target_len[i]));
    RETVAL = (SV*) newRV((SV*)avref); 
  OUTPUT:
    RETVAL

SV*
bam_text(bamh, ...)
  Bio::DB::Bam::Header bamh
  PREINIT:
    char   *newtext;
    STRLEN n;
  CODE:
    /* in case text is not null terminated, we copy it */
    RETVAL = newSVpv(bamh->text,bamh->l_text);
    if (items > 1) {
      newtext = (char*) SvPV(ST(1),n);
      strcpy(bamh->text,newtext);
      bamh->l_text = n;
    }
  OUTPUT:
    RETVAL

void
bam_parse_region(bamh,region)
    Bio::DB::Bam::Header bamh
    char*            region
    PROTOTYPE: $
    PREINIT:
       int seqid,start,end;
    PPCODE:
    {
      bam_parse_region(bamh,
		       region,
		       &seqid,
		       &start,
		       &end);
      if (seqid < 0)
	XSRETURN_EMPTY;
      else {
	EXTEND(sp,3);
	PUSHs(sv_2mortal(newSViv(seqid)));
	PUSHs(sv_2mortal(newSViv(start)));
	PUSHs(sv_2mortal(newSViv(end)));
      }
    }

void
bam_view1(bamh,alignment)
     Bio::DB::Bam::Header     bamh
     Bio::DB::Bam::Alignment  alignment
     PROTOTYPE: $$
     CODE:
       bam_view1(bamh,alignment);

void
bam_DESTROY(bamh)
  Bio::DB::Bam::Header bamh
  PROTOTYPE: $
  CODE:
    bam_header_destroy(bamh);

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Bam::Index PREFIX=bami_

int
bami_fetch(bai,bfp,ref,start,end,callback,callbackdata=&PL_sv_undef)
  Bio::DB::Bam::Index bai
  Bio::DB::Bam        bfp
  int   ref
  int   start
  int   end
  CV*   callback
  SV*   callbackdata
PREINIT:
  fetch_callback_data fcd;
CODE:
  {
    fcd.callback = (SV*) callback;
    fcd.data     = callbackdata;
    RETVAL = bam_fetch(bfp,bai,ref,start,end,&fcd,bam_fetch_fun);
  }
OUTPUT:
    RETVAL

void
bami_lpileup(bai,bfp,ref,start,end,callback,callbackdata=&PL_sv_undef)
  Bio::DB::Bam::Index bai
  Bio::DB::Bam        bfp
  int   ref
  int   start
  int   end
  CV*   callback
  SV*   callbackdata
PREINIT:  
  fetch_callback_data fcd;
  bam_lplbuf_t        *pileup;
CODE:
  fcd.callback = (SV*) callback;
  fcd.data     = callbackdata;
  pileup       = bam_lplbuf_init(invoke_pileup_callback_fun,(void*)&fcd);
  bam_fetch(bfp,bai,ref,start,end,(void*)pileup,add_lpileup_line);
  bam_lplbuf_push(NULL,pileup);
  bam_lplbuf_destroy(pileup);

void
bami_pileup(bai,bfp,ref,start,end,callback,callbackdata=&PL_sv_undef)
  Bio::DB::Bam::Index bai
  Bio::DB::Bam        bfp
  int   ref
  int   start
  int   end
  CV*   callback
  SV*   callbackdata
PREINIT:  
  fetch_callback_data fcd;
  bam_plbuf_t        *pileup;
CODE:
  fcd.callback = (SV*) callback;
  fcd.data     = callbackdata;
  pileup       = bam_plbuf_init(invoke_pileup_callback_fun,(void*)&fcd);
  bam_plp_set_maxcnt(pileup->iter,MaxPileupCnt);	
  bam_fetch(bfp,bai,ref,start,end,(void*)pileup,add_pileup_line);
  bam_plbuf_push(NULL,pileup);
  bam_plbuf_destroy(pileup);

AV*
bami_coverage(bai,bfp,ref,start,end,bins=0,maxcnt=8000)
    Bio::DB::Bam::Index bai
    Bio::DB::Bam        bfp
    int             ref
    int             start
    int             end
    int             bins
    int             maxcnt
PREINIT:
    coverage_graph  cg;
    bam_plbuf_t    *pileup;
    AV*             array;
    SV*             cov;
    int             i;
    bam_header_t   *bh;
CODE:
  {
 
      if (end >= MAX_REGION) {
          bgzf_seek(bfp,0,0);
          bh  = bam_header_read(bfp);
          end = bh->target_len[ref];
          bam_header_destroy(bh);
      }
      if ((bins==0) || (bins > (end-start)))
         bins = end-start;

      /* coverage graph used to communicate to our callback
	  the region we are sampling */
      cg.start = start;
      cg.end   = end;
      cg.reads = 0;
      cg.width = ((double)(end-start))/bins;
      Newxz(cg.bin,bins+1,int);

      /* accumulate coverage into the coverage graph */
      pileup   = bam_plbuf_init(coverage_from_pileup_fun,(void*)&cg);
      if (items >= 7)
            bam_plp_set_maxcnt(pileup->iter,maxcnt);
      else
            bam_plp_set_maxcnt(pileup->iter,MaxPileupCnt);
      bam_fetch(bfp,bai,ref,start,end,(void*)pileup,add_pileup_line);
      bam_plbuf_push(NULL,pileup);
      bam_plbuf_destroy(pileup);

      /* now normalize to coverage/bp and convert into an array */
      array = newAV();
      av_extend(array,bins);
      for  (i=0;i<bins;i++)
           av_store(array,i,newSVnv(((float)cg.bin[i])/cg.width));
      Safefree(cg.bin);
      RETVAL = array;
      sv_2mortal((SV*)RETVAL);  /* this fixes a documented bug in perl typemap */

  }
OUTPUT:
    RETVAL

void
bami_DESTROY(bai)
  Bio::DB::Bam::Index bai
  CODE:
    bam_index_destroy(bai);

MODULE = Bio::DB::Sam PACKAGE = Bio::DB::Bam::Pileup PREFIX=pl_

int
pl_qpos(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->qpos;
  OUTPUT:
    RETVAL

int
pl_pos(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->qpos+1;
  OUTPUT:
    RETVAL

int
pl_indel(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->indel;
  OUTPUT:
    RETVAL

int
pl_level(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->level;
  OUTPUT:
    RETVAL

int
pl_is_del(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->is_del;
  OUTPUT:
    RETVAL

int
pl_is_refskip(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->is_refskip;
  OUTPUT:
    RETVAL

int
pl_is_head(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->is_head;
  OUTPUT:
    RETVAL

int
pl_is_tail(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = pl->is_tail;
  OUTPUT:
    RETVAL

Bio::DB::Bam::Alignment
pl_b(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = bam_dup1(pl->b);
  OUTPUT:
     RETVAL

Bio::DB::Bam::Alignment
pl_alignment(pl)
  Bio::DB::Bam::Pileup pl
  CODE:
    RETVAL = bam_dup1(pl->b);
  OUTPUT:
     RETVAL