#include <stdlib.h>
 #include <string.h>
 
 #include <gstruct/bool.h>
 #include <gstruct/datadir.h>
 #include <gstruct/iit-read.h>
 #include <gstruct/genome.h>
 #include <gstruct/interval.h>
 #include <gstruct/genomicpos.h>
 #include <gstruct/complement.h>
 
 #include "gmapR.h"
 #include "genome.h"
 
 IIT_T
 readChromosomeIIT (const char *genome_dir, const char *db) {
   char *genomesubdir, *fileroot, *dbversion, *iitfile;
   genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,
                                            (char *)genome_dir, (char *)db);
   
   iitfile = (char *) calloc(strlen(genomesubdir) + strlen("/") +
                             strlen(fileroot) + strlen(".chromosome.iit") + 1,
                             sizeof(char));
   sprintf(iitfile, "%s/%s.chromosome.iit", genomesubdir, fileroot);
   IIT_T chromosome_iit = IIT_read(iitfile, /*name*/NULL, /*readonlyp*/true,
                                   /*divread*/READ_ALL, /*divstring*/NULL,
                                   /*add_iit_p*/false, /*labels_read_p*/true);
   free(iitfile);
   free(fileroot);
   free(dbversion);
   free(genomesubdir);
 
   return chromosome_iit;
 }
 
 Genome_T
 createGenome (const char *genome_dir, const char *db) {
   char *genomesubdir, *fileroot, *dbversion;
   genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,
                                            (char *)genome_dir, (char *)db);
   
   Genome_T genome = Genome_new(genomesubdir, fileroot, /*snps_root*/NULL,
                                /*uncompressedp*/false, /*access*/USE_MMAP_ONLY);
   free(fileroot);
   free(dbversion);
   free(genomesubdir);
 
   return genome;
 }
 
 static char complCode[128] = COMPLEMENT_LC;
 
 static void
 make_complement_inplace (char *sequence, Genomicpos_T length) {
   char temp;
   int i, j;
 
   for (i = 0, j = length-1; i < length/2; i++, j--) {
     temp = complCode[(int) sequence[i]];
     sequence[i] = complCode[(int) sequence[j]];
     sequence[j] = temp;
   }
   if (i == j) {
     sequence[i] = complCode[(int) sequence[i]];
   }
 
   return;
 }
 
 SEXP
 R_Genome_getSeq (SEXP genome_dir_R, SEXP db_R,
                  SEXP seqnames_R, SEXP start_R, SEXP width_R, SEXP strand_R)
 {
   const char *genome_dir =
     genome_dir_R == R_NilValue ? NULL : CHAR(asChar(genome_dir_R));
   const char *db = CHAR(asChar(db_R));
   Genome_T genome = createGenome(genome_dir, db);
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
   int *start = INTEGER(start_R);
   int *width = INTEGER(width_R);
 
   int max_width = 0;
   for (int i = 0; i < length(width_R); i++) {
     if (width[i] > max_width)
       max_width = width[i];
   }
   char *buffer = (char *) R_alloc(sizeof(char), max_width + 1);
 
   int seqname_index;
   const char *old_seqname = NULL; // premature optimization
   SEXP result;
   PROTECT(result = allocVector(STRSXP, length(start_R)));
   for (int i = 0; i < length(start_R); i++) {
     const char *seqname = CHAR(STRING_ELT(seqnames_R, i));
     if (old_seqname == NULL || strcmp(seqname, old_seqname)) {
       seqname_index = IIT_find_linear(chromosome_iit, (char *)seqname);
       if (seqname_index < 0) {
         error("Cannot find chromosome %s in genome", seqname);
       }
     }
     old_seqname = seqname;
     int chroffset = Interval_low(IIT_interval(chromosome_iit, seqname_index));
     Genome_fill_buffer_simple(genome, chroffset + start[i] - 1, width[i],
                               buffer);
     if (CHAR(STRING_ELT(strand_R, i))[0] == '-')
       make_complement_inplace(buffer, width[i]);
     SET_STRING_ELT(result, i, mkChar(buffer));
   }
 
   IIT_free(&chromosome_iit);
   Genome_free(&genome);
 
   UNPROTECT(1);
   return result;
 }