#define DEBUG2 0

#include <stdlib.h> /* for abs */
#include <string.h> /* for strlen */

#include "iit.h"
#include "bytestream.h"

enum { SEQNAMES, POS, REF, READ, N_CYCLES, N_CYCLES_REF, COUNT, COUNT_REF,
       RAW_COUNT_TOTAL, COUNT_TOTAL, COUNT_PLUS, COUNT_PLUS_REF, COUNT_MINUS,
       COUNT_MINUS_REF, DELCOUNT_PLUS, DELCOUNT_MINUS,
       READ_POS_MEAN, READ_POS_MEAN_REF, READ_POS_VAR,
       READ_POS_VAR_REF, MDFNE, MDFNE_REF,  CODON_STRAND,
       COUNT_XS_PLUS, COUNT_XS_PLUS_REF,
       COUNT_XS_MINUS, COUNT_XS_MINUS_REF,
       COUNT_HIGH_NM, COUNT_HIGH_NM_REF, N_BASE_COLS };

enum { STRAND_PLUS = 1,
       STRAND_MINUS,
       STRAND_NONE };

static char *codon_table[64] = 
  {"AAA", "AAC", "AAG", "AAT", "ACA", "ACC", "ACG", "ACT",
   "AGA", "AGC", "AGG", "AGT", "ATA", "ATC", "ATG", "ATT",
   "CAA", "CAC", "CAG", "CAT", "CCA", "CCC", "CCG", "CCT",
   "CGA", "CGC", "CGG", "CGT", "CTA", "CTC", "CTG", "CTT",
   "GAA", "GAC", "GAG", "GAT", "GCA", "GCC", "GCG", "GCT",
   "GGA", "GGC", "GGG", "GGT", "GTA", "GTC", "GTG", "GTT",
   "TAA", "TAC", "TAG", "TAT", "TCA", "TCC", "TCG", "TCT",
   "TGA", "TGC", "TGG", "TGT", "TTA", "TTC", "TTG", "TTT"};

static int NUM_PTRS = 5;

typedef struct TallyTable {
  SEXP seqnames_R;
  int *pos;
  SEXP ref_R;
  SEXP read_R;
  int *n_cycles;
  int *n_cycles_ref;
  int *count;
  int *count_ref;
  int *raw_count_total;
  int *count_total;
  int *count_plus;
  int *count_plus_ref;
  int *count_minus;
  int *count_minus_ref;
  int *delcount_plus;
  int *delcount_minus;
  double *read_pos_mean;
  double *read_pos_mean_ref;
  double *read_pos_var;
  double *read_pos_var_ref;
  double *mdfne;
  double *mdfne_ref;
  int *count_xs_plus;
  int *count_xs_plus_ref;
  int *count_xs_minus;
  int *count_xs_minus_ref;
  int *count_high_nm;
  int *count_high_nm_ref;
  int *strand;
  int **cycle_bins;
} TallyTable;

typedef struct TallyParam {
  int *cycle_breaks;
  int n_cycle_bins;
  int read_length;
  double *mdfne_buf;
  bool xs;
  int high_nm_score;
} TallyParam;

static SEXP R_TallyTable_new(int n_rows, int n_cycle_bins, bool xs) {
  SEXP tally_R; /* the result list */
  PROTECT(tally_R = allocVector(VECSXP, N_BASE_COLS + n_cycle_bins));
  SET_VECTOR_ELT(tally_R, SEQNAMES, allocVector(STRSXP, n_rows));
  SET_VECTOR_ELT(tally_R, POS, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, REF, allocVector(STRSXP, n_rows));
  SET_VECTOR_ELT(tally_R, READ, allocVector(STRSXP, n_rows));
  SET_VECTOR_ELT(tally_R, N_CYCLES, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, N_CYCLES_REF, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_REF, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, RAW_COUNT_TOTAL, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_TOTAL, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_PLUS, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_PLUS_REF, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_MINUS, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_MINUS_REF, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, READ_POS_MEAN, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, READ_POS_MEAN_REF, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, READ_POS_VAR, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, READ_POS_VAR_REF, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, MDFNE, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, MDFNE_REF, allocVector(REALSXP, n_rows));
  SET_VECTOR_ELT(tally_R, CODON_STRAND, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, DELCOUNT_PLUS, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, DELCOUNT_MINUS, allocVector(INTSXP, n_rows));
  if (xs) {
      SET_VECTOR_ELT(tally_R, COUNT_XS_PLUS, allocVector(INTSXP, n_rows));
      SET_VECTOR_ELT(tally_R, COUNT_XS_PLUS_REF, allocVector(INTSXP, n_rows));
      SET_VECTOR_ELT(tally_R, COUNT_XS_MINUS, allocVector(INTSXP, n_rows));
      SET_VECTOR_ELT(tally_R, COUNT_XS_MINUS_REF, allocVector(INTSXP, n_rows));
  }
  SET_VECTOR_ELT(tally_R, COUNT_HIGH_NM, allocVector(INTSXP, n_rows));
  SET_VECTOR_ELT(tally_R, COUNT_HIGH_NM_REF, allocVector(INTSXP, n_rows));
  
  for (int bin = 0; bin < n_cycle_bins; bin++) {
    SEXP cycle_bin_R = allocVector(INTSXP, n_rows);
    SET_VECTOR_ELT(tally_R, bin + N_BASE_COLS, cycle_bin_R);
  }

  UNPROTECT(1);
  return tally_R;
}

static TallyTable *TallyTable_new(SEXP tally_R, bool xs) {
  TallyTable *tally = (TallyTable *) S_alloc(sizeof(TallyTable), 1);
  int n_cycle_bins = length(tally_R) - N_BASE_COLS;
  
  tally->seqnames_R = VECTOR_ELT(tally_R, SEQNAMES);
  tally->pos = INTEGER(VECTOR_ELT(tally_R, POS));
  tally->ref_R = VECTOR_ELT(tally_R, REF);
  tally->read_R = VECTOR_ELT(tally_R, READ);
  tally->n_cycles = INTEGER(VECTOR_ELT(tally_R, N_CYCLES));
  tally->n_cycles_ref = INTEGER(VECTOR_ELT(tally_R, N_CYCLES_REF));
  tally->count = INTEGER(VECTOR_ELT(tally_R, COUNT));
  tally->count_ref = INTEGER(VECTOR_ELT(tally_R, COUNT_REF));
  tally->raw_count_total = INTEGER(VECTOR_ELT(tally_R, RAW_COUNT_TOTAL));
  tally->count_total = INTEGER(VECTOR_ELT(tally_R, COUNT_TOTAL));
  tally->count_plus = INTEGER(VECTOR_ELT(tally_R, COUNT_PLUS));
  tally->count_plus_ref = INTEGER(VECTOR_ELT(tally_R, COUNT_PLUS_REF));
  tally->count_minus = INTEGER(VECTOR_ELT(tally_R, COUNT_MINUS));
  tally->count_minus_ref = INTEGER(VECTOR_ELT(tally_R, COUNT_MINUS_REF));
  tally->read_pos_mean = REAL(VECTOR_ELT(tally_R, READ_POS_MEAN));
  tally->read_pos_mean_ref = REAL(VECTOR_ELT(tally_R, READ_POS_MEAN_REF));
  tally->read_pos_var = REAL(VECTOR_ELT(tally_R, READ_POS_VAR));
  tally->read_pos_var_ref = REAL(VECTOR_ELT(tally_R, READ_POS_VAR_REF));
  tally->mdfne = REAL(VECTOR_ELT(tally_R, MDFNE));
  tally->mdfne_ref = REAL(VECTOR_ELT(tally_R, MDFNE_REF));
  tally->strand = INTEGER(VECTOR_ELT(tally_R, CODON_STRAND));
  tally->delcount_plus = INTEGER(VECTOR_ELT(tally_R, DELCOUNT_PLUS));
  tally->delcount_minus = INTEGER(VECTOR_ELT(tally_R, DELCOUNT_MINUS));
  tally->cycle_bins = (int **) R_alloc(sizeof(int*), n_cycle_bins);
  if (xs) {
      tally->count_xs_plus = INTEGER(VECTOR_ELT(tally_R, COUNT_XS_PLUS));
      tally->count_xs_plus_ref = INTEGER(VECTOR_ELT(tally_R,
                                                    COUNT_XS_PLUS_REF));
      tally->count_xs_minus = INTEGER(VECTOR_ELT(tally_R, COUNT_XS_MINUS));
      tally->count_xs_minus_ref = INTEGER(VECTOR_ELT(tally_R,
                                                     COUNT_XS_MINUS_REF));
  }
  tally->count_high_nm = INTEGER(VECTOR_ELT(tally_R, COUNT_HIGH_NM));
  tally->count_high_nm_ref = INTEGER(VECTOR_ELT(tally_R,
                                                COUNT_HIGH_NM_REF));
  
  for (int bin = 0; bin < n_cycle_bins; bin++) {
    tally->cycle_bins[bin] = INTEGER(VECTOR_ELT(tally_R, bin + N_BASE_COLS));
  }

  return tally;
}

static void
read_total_counts(unsigned char **bytes, int row, int *raw_count_total,
                  int *count_total)
{
  raw_count_total[row] = read_int(bytes);
  int count_total_plus = read_int(bytes);
  int count_total_minus = read_int(bytes);
  count_total[row] = count_total_plus + count_total_minus;
}

static void
read_nm_counts(unsigned char **bytes, int row, int *count_high_nm,
               int high_nm_score)
{
    int n_nm = read_int(bytes);
    count_high_nm[row] = 0;
    /* FIXME: currently bam_tally gives us counts whether we want them or not */
    for (int index = 0; index < n_nm; index++) {
        int nm = read_int(bytes);
	int count = read_int(bytes);
        if (nm >= high_nm_score) {
            count_high_nm[row] += count;
        }
    }
}

static void
read_xs_counts(unsigned char **bytes, int row, int *count_xs_plus,
               int *count_xs_minus)
{
    int n_xs = read_int(bytes);
    if (count_xs_plus == NULL) {
        return;
    }
    count_xs_plus[row] = 0;
    count_xs_minus[row] = 0;
    for (int index = 0; index < n_xs; index++) {
	int xs = read_int(bytes);
	int count = read_int(bytes);
	if (xs == 1) {
	    count_xs_plus[row] = count;
	} else if (xs == 2) {
	    count_xs_minus[row] = count;
	} 
    }
}

static void
read_cycle_counts(unsigned char **bytes, int row, TallyParam param,
                  int *n_cycles, double *read_pos_mean, double *read_pos_var,
                  double *mdfne, int **cycle_bins)
{
  int n_cycle_breaks = param.n_cycle_bins + 1;
  int count_sum = 0, weighted_sum = 0, weighted_sum_sq = 0;
  int midpoint = param.read_length / 2.0 + 0.5;

  if (param.mdfne_buf != NULL) {
    memset(param.mdfne_buf, 0, sizeof(int) * midpoint);
  }
  n_cycles[row] = read_int(bytes);
  
  for (int index = 0; index < n_cycles[row]; index++) {
    int cycle = abs(read_int(bytes));
    int count = read_int(bytes);
    int bin = 0;
    if (param.mdfne_buf != NULL && cycle <= param.read_length) {
      param.mdfne_buf[midpoint-abs(cycle-midpoint)-1] = count;
    }
    count_sum += count;
    weighted_sum += cycle * count;
    weighted_sum_sq += cycle * cycle * count;
    if (param.n_cycle_bins > 0) {
      while(n_cycle_breaks > bin &&
            cycle > param.cycle_breaks[bin])
        bin++;
      if (bin > 0 && bin < n_cycle_breaks) {
        cycle_bins[bin-1][row] += count;
      }
    }
  }

  read_pos_mean[row] = ((double)weighted_sum) / count_sum;
  if (count_sum > 1) {
    read_pos_var[row] = ((double)weighted_sum_sq) / (count_sum - 1) -
      (count_sum / (count_sum - 1)) * read_pos_mean[row] * read_pos_mean[row];
  } else {
    read_pos_var[row] = NA_REAL;
  }
  
  mdfne[row] = NA_REAL;
  if (param.mdfne_buf != NULL) {
    int prev_dist = 0;
    int total_count = count_sum;
    count_sum = 0;
    for (int dist = 0; dist < midpoint; dist++) {
      count_sum += param.mdfne_buf[dist];
      if (count_sum > total_count / 2) {
        if (total_count % 2 == 0) {
          mdfne[row] = prev_dist + (dist - prev_dist) / 2.0;
        } else {
          mdfne[row] = dist;
        }
        break;
      }
      if (param.mdfne_buf[dist] > 0) {
        prev_dist = dist;
      }
    }
  }
}

static int
parse_indels(unsigned char *bytes, int row,
             TallyParam param, TallyTable *tally, bool insertion)
{
  int indel_count = read_int(&bytes);
  for (int indel = 0; indel < indel_count; indel++, row++) {
    for (int b = 0; b < param.n_cycle_bins; b++) {
      tally->cycle_bins[b][row] = 0;
    }
    tally->count_plus[row] = read_int(&bytes);
    tally->count_minus[row] = read_int(&bytes);
    tally->count[row] = tally->count_plus[row] + tally->count_minus[row];
    tally->count_plus_ref[row] = read_int(&bytes);
    tally->count_minus_ref[row] = read_int(&bytes);
    tally->count_ref[row] = tally->count_plus_ref[row] +
	tally->count_minus_ref[row];
    tally->raw_count_total[row] = tally->count_ref[row] + tally->count[row];
    tally->count_total[row] = tally->raw_count_total[row];
    tally->strand[row] = STRAND_NONE;
    SEXP seq_R = mkChar(read_string(&bytes));
    if (insertion) {
      SET_STRING_ELT(tally->read_R, row, seq_R);
      SET_STRING_ELT(tally->ref_R, row, R_BlankString);
    }
    else {
      SET_STRING_ELT(tally->read_R, row, R_BlankString);
      SET_STRING_ELT(tally->ref_R, row, seq_R);
    }
    read_cycle_counts(&bytes, row, param, tally->n_cycles,
                      tally->read_pos_mean, tally->read_pos_var,
                      tally->mdfne, tally->cycle_bins);
    read_nm_counts(&bytes, row, tally->count_high_nm, param.high_nm_score);
    read_xs_counts(&bytes, row, tally->count_xs_plus, tally->count_xs_minus);
    /* no position from which to tabulate the cycles */
    tally->n_cycles_ref[row] = NA_INTEGER;
    tally->read_pos_mean_ref[row] = NA_REAL;
    tally->read_pos_var_ref[row] = NA_REAL;
    tally->mdfne_ref[row] = NA_REAL;
    /* overlapping deletion count only relevant for SNVs */
    tally->delcount_plus[row] = NA_INTEGER;
    tally->delcount_minus[row] = NA_INTEGER;
    /* information not available for "ref" allele */
    tally->count_high_nm_ref[row] = NA_INTEGER;
    if (param.xs) {
        tally->count_xs_plus_ref[row] = NA_INTEGER;
        tally->count_xs_minus_ref[row] = NA_INTEGER;
    }
  }
  return indel_count;
}

static int
read_allele_counts(unsigned char **bytes, int row, SEXP read_R,
                   int *count_plus, int *count_minus, int *count,
                   int *delcount_plus, int *delcount_minus,
                   int strand)
{
  int n_alleles = 0;
  char allele;
  char stop = strand == STRAND_NONE ? '\0' : (char) 255;
  
  while((allele = (char)read_char(bytes)) != stop) {
    if (allele == '_') {
        *delcount_plus = read_int(bytes);
        *delcount_minus = read_int(bytes);
        continue;
    }
    if(strand == STRAND_NONE)
       SET_STRING_ELT(read_R, row, mkCharLen(&allele, 1));
    else
       SET_STRING_ELT(read_R, row, mkCharLen(codon_table[(int)allele], 3));
    count_plus[row] = read_int(bytes);
    count_minus[row] = read_int(bytes);
    count[row] = count_plus[row] + count_minus[row];
    row++;
    n_alleles++;
  }

  return n_alleles;
}

static int
parse_alleles(unsigned char *bytes, int row, int ref_row,
              TallyParam param, TallyTable *tally, int strand)
{
  read_total_counts(&bytes, row, tally->raw_count_total, tally->count_total);
  int delcount_plus = 0, delcount_minus = 0;
  int n_alleles = read_allele_counts(&bytes, row, tally->read_R,
                                     tally->count_plus, tally->count_minus,
	                             tally->count,
                                     &delcount_plus, &delcount_minus,
	                             strand);
  for (int allele = 0; allele < n_alleles; allele++, row++) {
    tally->n_cycles[row] = 0;
    for (int b = 0; b < param.n_cycle_bins; b++) {
      tally->cycle_bins[b][row] = 0;
    }
    tally->strand[row] = strand;
    tally->delcount_plus[row] = delcount_plus;
    tally->delcount_minus[row] = delcount_minus;
    if (tally->count[row] > 0) {
      read_cycle_counts(&bytes, row, param, tally->n_cycles,
                        tally->read_pos_mean, tally->read_pos_var,
                        tally->mdfne, tally->cycle_bins);
      read_nm_counts(&bytes, row, tally->count_high_nm, param.high_nm_score);
      if (param.xs) {
        read_xs_counts(&bytes, row, tally->count_xs_plus,
                       tally->count_xs_minus);
      }
    } else {
      tally->read_pos_mean[row] = R_NaN;
      tally->read_pos_var[row] = NA_REAL;
      tally->mdfne[row] = NA_REAL;
      tally->count_high_nm[row] = 0;
      if (param.xs) {
        tally->count_xs_plus[row] = 0;
        tally->count_xs_minus[row] = 0;
      }
    }
  }
  for (int r = ref_row; r < row; r++) {
    tally->n_cycles_ref[r] = tally->n_cycles[ref_row];
    tally->raw_count_total[r] = tally->raw_count_total[ref_row];
    tally->count_total[r] = tally->count_total[ref_row];
    tally->count_plus_ref[r] = tally->count_plus[ref_row];
    tally->count_minus_ref[r] = tally->count_minus[ref_row];
    tally->count_ref[r] = tally->count_plus_ref[r] + tally->count_minus_ref[r];
    tally->read_pos_mean_ref[r] = tally->read_pos_mean[ref_row];
    tally->read_pos_var_ref[r] = tally->read_pos_var[ref_row];
    tally->mdfne_ref[r] = tally->mdfne[ref_row];
    tally->count_high_nm_ref[r] = tally->count_high_nm[ref_row];
    if (param.xs) {
      tally->count_xs_plus_ref[r] = tally->count_xs_plus[ref_row];
      tally->count_xs_minus_ref[r] = tally->count_xs_minus[ref_row];
    }
    SET_STRING_ELT(tally->ref_R, r, STRING_ELT(tally->read_R, ref_row));  
  }

    /* clear the 'alt' columns for the 'ref' row with NAs */
    SET_STRING_ELT(tally->read_R, ref_row, NA_STRING);
    tally->n_cycles[ref_row] = NA_INTEGER;
    tally->count_plus[ref_row] = NA_INTEGER;
    tally->count_minus[ref_row] = NA_INTEGER;
    tally->count[ref_row] = NA_INTEGER;
    tally->read_pos_mean[ref_row] = NA_REAL;
    tally->read_pos_var[ref_row] = NA_REAL;
    tally->mdfne[ref_row] = NA_REAL;
    if (param.xs) {
      tally->count_xs_plus[ref_row] = NA_INTEGER;
      tally->count_xs_minus[ref_row] = NA_INTEGER;
    }
    
  return n_alleles;
}
    
static int parse_indel_count(unsigned char *bytes) {
  int count = read_int(&bytes);
  return count;
}

static int parse_allele_count(unsigned char *bytes) {
  int n_alleles = 1; /* always have a reference */
  bytes += sizeof(int) * 5 + 1; /* skip total and reference */
  while(bytes[0] != '\0') {
      if (bytes[0] != '_')
          n_alleles++;
      bytes += sizeof(int) * 2 + 1;
  }
  return n_alleles;
}


static int parse_codon_count(unsigned char *bytes) {
  int n_alleles = 1; /* always have a reference */
  bytes += sizeof(int) * 4 + 1; /* skip reference */
  while((int) bytes[0] != 255 ){ // && bytes [0] != '\0') {
    bytes += sizeof(int) * 2 + 1;
    n_alleles++;
  }
  return n_alleles;
}


static int count_rows_for_interval(IIT_T tally_iit, int index) {
  int n_rows = 0;
  unsigned char *bytes = IIT_data(tally_iit, index);
  int width = IIT_length(tally_iit, index);
  unsigned char *base = bytes + (NUM_PTRS * width + 1) * sizeof(int);
  for (int pos = 0; pos < width; pos++) {
    int insertion_offset = read_int(&bytes);
    int deletion_offset = read_int(&bytes);
    int allele_offset = read_int(&bytes);
    int codon_plus_offset = read_int(&bytes);
    int codon_minus_offset = read_int(&bytes);
    int next_offset = read_int(&bytes);
    bytes -= 4; /* rewind from read-ahead */
    if (deletion_offset - insertion_offset > 0) {
      n_rows += parse_indel_count(base + insertion_offset);
    }
    if (allele_offset - deletion_offset > 0) {
      n_rows += parse_indel_count(base + deletion_offset);
    }
    if (codon_plus_offset - allele_offset > 0) {
      n_rows += parse_allele_count(base + allele_offset);
    }
    if (codon_minus_offset - codon_plus_offset > 0) {
      n_rows += parse_codon_count(base + codon_plus_offset);
    }
    if (next_offset - codon_minus_offset > 0) {
      n_rows += parse_codon_count(base + codon_minus_offset);
    } 

  }
  return n_rows;
}

static int parse_interval(IIT_T tally_iit, int index,
                          TallyParam param, TallyTable *tally, int row)
{
  unsigned char *bytes = IIT_data(tally_iit, index);
  int width = IIT_length(tally_iit, index);
  unsigned char *base = bytes + (NUM_PTRS * width + 1) * sizeof(int);
  int start = IIT_interval_low(tally_iit, index);
  SEXP divstring_R;
  PROTECT(divstring_R = mkChar(IIT_divstring_from_index(tally_iit, index)));
  for (int position = 0; position < width; position++) {
    int insertion_offset = read_int(&bytes);
    int deletion_offset = read_int(&bytes);
    int allele_offset = read_int(&bytes);
    int codon_plus_offset = read_int(&bytes);
    int codon_minus_offset = read_int(&bytes);
    int next_offset = read_int(&bytes);
    int ref_row = row;
    bytes -= 4; /* rewind from read-ahead */
    if (codon_plus_offset - allele_offset > 0)
      row += parse_alleles(base + allele_offset, row, ref_row,
                           param, tally, STRAND_NONE);
    if (deletion_offset - insertion_offset > 0) 
      row += parse_indels(base + insertion_offset, row,
                          param, tally, true);
    if (allele_offset - deletion_offset > 0)
      row += parse_indels(base + deletion_offset, row,
                          param, tally, false);

    if (codon_minus_offset - codon_plus_offset > 0) {
      row += parse_alleles(base + codon_plus_offset, row, row,
                           param, tally, STRAND_PLUS);
    }
    if (next_offset - codon_minus_offset > 0) { 
      row += parse_alleles(base + codon_minus_offset, row, row,
                           param, tally, STRAND_MINUS);
    }
    /* fill in position information */
    for (int r = ref_row; r < row; r++) {
      SET_STRING_ELT(tally->seqnames_R, r, divstring_R);
      tally->pos[r] = start + position;
    }      
  }
  UNPROTECT(1);
  return row;
}

static SEXP parse_all(IIT_T tally_iit, TallyParam param)
{
  int n_rows = 0;
  /* loop over the IIT, getting the total number of rows
     this is (num alts + 1) for every position. */
  for (int index = 1; index <= IIT_total_nintervals(tally_iit); index++) {
    n_rows += count_rows_for_interval(tally_iit, index);
  }
  
  SEXP tally_R;
  PROTECT(tally_R = R_TallyTable_new(n_rows, param.n_cycle_bins, param.xs));
  TallyTable *tally = TallyTable_new(tally_R, param.xs);
  
  int row = 0;
  for (int index = 1; index <= IIT_total_nintervals(tally_iit); index++) {
    row = parse_interval(tally_iit, index, param, tally, row);
  }

  UNPROTECT(1);
  return tally_R;
}

static SEXP parse_some(IIT_T tally_iit, TallyParam param,
                       SEXP chr_R, int *start, int *end)
{
  int n_rows = 0;
  int *nmatches = (int *) R_alloc(sizeof(int), length(chr_R));
  int **indexes = (int **) R_alloc(sizeof(int*), length(chr_R));
  
  for (int i = 0; i < length(chr_R); i++) {
    indexes[i] = IIT_get(nmatches + i, tally_iit,
                         (char *) CHAR(STRING_ELT(chr_R, i)),
                         start[i], end[i], false);
    for (int j = 0; j < nmatches[i]; j++) {
      n_rows += count_rows_for_interval(tally_iit, indexes[i][j]);
    }
  }
  
  SEXP tally_R;
  PROTECT(tally_R = R_TallyTable_new(n_rows, param.n_cycle_bins, param.xs));
  TallyTable *tally = TallyTable_new(tally_R, param.xs);

  int row = 0;
  for (int i = 0; i < length(chr_R); i++) {
    for (int j = 0; j < nmatches[i]; j++) {
      row = parse_interval(tally_iit, indexes[i][j], param, tally, row);
    }
  }
  
  UNPROTECT(1);
  return tally_R;
}

/* FORMAT

   block header:
   [0][ins][del][mm][codon+][codon-]

   mismatches/codons:
   [t+][t-][rn][r+][r-]([an][a+][a-])*
   mismatches only:
   '_'[d+][d-][\0]([c#]([cv][cc])*[nm#]([nmv][nmc])*[xs#]([xsv][xsc])*)*

   insertions:
   [i#]([t+][t-][r+][r-][seq][c#]([cv][cc])*[nm#]([nmv][nmc])*[xs#]([xsv][xsc])*

   deletions:
   [d#]([t+][t-][r+][r-][seq][c#]([cv][cc])*[nm#]([nmv][nmc])*[xs#]([xsv][xsc])*
*/

SEXP R_tally_iit_parse(SEXP tally_iit_R, SEXP cycle_breaks_R,
                       SEXP which_R,
                       SEXP read_length_R, SEXP xs_R, SEXP high_nm_score_R)
{
  IIT_T tally_iit = (IIT_T) R_ExternalPtrAddr(tally_iit_R);
  SEXP tally_R;
  TallyParam param;
  param.cycle_breaks =
    cycle_breaks_R == R_NilValue ? NULL : INTEGER(cycle_breaks_R);
  param.n_cycle_bins =
    length(cycle_breaks_R) == 0 ? 0 : length(cycle_breaks_R) - 1;
  param.read_length = asInteger(read_length_R);
  if (param.read_length != NA_INTEGER) {
    param.mdfne_buf = (double *)R_alloc(sizeof(double), param.read_length);
  } else {
    param.mdfne_buf = NULL;
  }
  param.xs = asLogical(xs_R);
  param.high_nm_score = asInteger(high_nm_score_R);
  
  if (which_R == R_NilValue) {
    tally_R = parse_all(tally_iit, param);
  } else {
    SEXP chr_R = VECTOR_ELT(which_R, 0);
    int *start = INTEGER(VECTOR_ELT(which_R, 1));
    int *end = INTEGER(VECTOR_ELT(which_R, 2));
    tally_R = parse_some(tally_iit, param, chr_R, start, end);
  }
  
  return tally_R;
}