/*
  chip_part.c

  Partitioning Tool.
  The program implements an algorithm for large-scale 
  Chip-Seq data sets. 
  It heuristically partitions the data by using a 
  two-way bisection of the ChIPseq count frequencies 
  or density over the genome. 
  
  # Arguments:
  # feature type, density threshold,
  # transition penalty, counts threshold

  Giovanna Ambrosini, EPFL, Giovanna.Ambrosini@epfl.ch

  Copyright (c) 2014 EPFL and Swiss Institute of Bioinformatics.

  This program 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.

  This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

*/
//#define DEBUG 
#define _GNU_SOURCE
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <ctype.h>
#include <sys/stat.h>
#ifdef DEBUG
#include <mcheck.h>
#endif

#include "version.h"

/*#define BUF_SIZE 4096 */
#define BUF_SIZE 8192
#define LINE_SIZE 1024
#define FT_MAX  64
#define SEQ_ID  32
#define FT_MAX  64
#define SEQ_ID  32
#define POS_MAX 16
#define CNT_MAX 16
#define EXT_MAX 256

typedef struct _options_t {
  int help;
  int debug;
} options_t;

static options_t options;

typedef struct _feature_t {
   char seq_id[SEQ_ID];
   char *ft;
   char ft_str;
   char **name;
   char *strand;
   int *pos;
   int *cnt;
   int *npo;
   int *nct;
} feature_t, *feature_p_t;

feature_t ref_ft;
int strand_flag = 0;

char *Feature = NULL;
float Dthres = 0;
int Tpen = 0;
int Coff = 1;

unsigned long TotLen = 0;
unsigned long TotCnts = 0;
int TotSeqs = 0;
int NumSegs = 0;
unsigned long TotFragLen = 0;
unsigned long AvFragLen = 0;
unsigned long AvCntsFrag = 0;
unsigned long TotCntsFrag = 0;


void 
split_seq(int len, int end, int last_seq)
{
  /* Partitioning Algorithm  */
  /* X(i) partitioning ends in over-representation state 
     Y(i) partitioning ends in under-representation state
     q(i) trace-back codes

     Recursion

      X(i) = max X(i-1) + cnt(i) - (pos(i)-pos(i-1))*Dthres,                     (XX)
                 Y(i-1) + Tpen + cnt(i) - Dthres + (pos(i)-pos(i-1)-1)*Dthres,   (YX)
		 X(i-1) + 2*Tpen + cnt(i) - Dthres + (pos(i)-pos(i-1)-2)*Dthres  (ZX)

      Y(i) = max X(i-1) + Tpen + (pos(i)-pos(i-1))*Dthres - cnt(i),              (XY)
                 Y(i-1) + (pos(i)-pos(i-1))*Dthres - cnt(i)                      (YY)
  */
  int i, j, k = 1;
  size_t mLen1 = BUF_SIZE;
  size_t mLen2 = BUF_SIZE;
  size_t mLen3 = BUF_SIZE;
  float Xmax = 0;
  float Ymax = 0;
  float XX, YX, ZX, YY, XY;
  int *q;
  int **Reg;
  int *cnt;
  int seq_len = 0;
  int tot_cnts = 0;

  if (( q = (int*)calloc(mLen1, sizeof(int))) == NULL) {
    perror("split_seq: malloc");
    exit(1);
  }
  if (( cnt = (int*)calloc(mLen2, sizeof(int))) == NULL) {
    perror("split_seq: malloc");
    exit(1);
  }
  if (( Reg = (int**)calloc((size_t)2, sizeof(int*))) == NULL) {
    perror("split_seq: malloc");
    exit(1);
  }
  for(i = 0; i < 2; i++) {
    Reg[i] = (int*)calloc(mLen3, sizeof(int));
  } 
  if (end == 0) { 
    ref_ft.pos[len + 1] = ref_ft.pos[len];
    ref_ft.name[len + 1] = ref_ft.name[len];
  } else 
    ref_ft.pos[len + 1] = end;

    ref_ft.cnt[len + 1] = 0;
    q[0] = 0;
  /* Fill in trace-back codes  q[i] */
  /* 
     0  Xmax = XX ; Ymax = YY
     1  Xmax = ZX ; Ymax = YY
     2  Xmax = YX ; Ymax = YY
     3  Xmax = XX ; Ymax = XY
     4  Xmax = ZX ; Ymax = XY
     5  Xmax = YX ; Ymax = XY
  */
  for (i = 1; i <= len + 1; i++) {
    XX = Xmax + ref_ft.cnt[i] - (ref_ft.pos[i] - ref_ft.pos[i - 1])*Dthres;
    YX = Ymax  + Tpen + ref_ft.cnt[i] - Dthres + (ref_ft.pos[i] - ref_ft.pos[i - 1] - 1)*Dthres;
    ZX = Xmax  + 2*Tpen + ref_ft.cnt[i] - Dthres + (ref_ft.pos[i] - ref_ft.pos[i - 1] - 1)*Dthres;
    YY = Ymax - ref_ft.cnt[i] + (ref_ft.pos[i] - ref_ft.pos[i - 1])*Dthres;
    XY = Xmax + Tpen - ref_ft.cnt[i] + (ref_ft.pos[i] - ref_ft.pos[i - 1])*Dthres;
    if ((unsigned int)i >= mLen1) {
      mLen1 *= 2;
      if (( q = (int *)realloc(q, mLen1 * sizeof(int))) == NULL) {
    	perror("split_seq: realloc");
    	exit(1);
      }
    }
    q[i] = 0;
    if ((YX < XX) && (ZX < XX)) { 
      Xmax = XX;
    } else if (YX < ZX) {
      Xmax = ZX;
      q[i] += 1;
    } else {
      Xmax = YX;
      q[i] += 2;
    }
    if (XY < YY) {
      Ymax = YY;
    } else {
      Ymax = XY;
      q[i] += 3;
    }
  }
  /* Trace-back */
  k = 1;
  if (Xmax > Ymax) {
    Reg[0][k] = 0;
    Reg[1][k] = ref_ft.pos[len + 1]; 
  } else {
    Reg[0][k] = 0;
    Reg[1][k] = 0; 
  }
  for (i = len + 1; i > 0; i--) {
    if (Reg[1][k] == 0) {
      if (q[i] > 2) 
        Reg[1][k] = ref_ft.pos[i - 1]; /* Enter a signal rich region 
					(we are in a poor region) */
    } else {
      if (q[i] > 2)
        q[i] -= 3;
      if ((unsigned int)k >= mLen2) {
        mLen2 *= 2;
        if ((cnt = (int*)realloc(cnt , mLen2 * sizeof(int))) == NULL) {
          perror("split_seq: realloc");
          exit(1);
        }
	memset((void *)&cnt[k], 0,  mLen2 * sizeof(int) / 2);
      }
      if (q[i] == 1) { /* Enter signal-rich region for the first time. 
      			We first close the previous one */
        Reg[0][k] = ref_ft.pos[i];
	cnt[k] += ref_ft.cnt[i];
	k++;
        if ((unsigned int)k >= mLen3) {
          mLen3 *= 2;
          for(j = 0; j < 2; j++) {
            if ((Reg[j] = (int*)realloc(Reg[j], mLen3 * sizeof(int))) == NULL) {
              perror("split_seq: realloc");
              exit(1);
	    }
          } 
        }
        Reg[0][k] = 0;
        Reg[1][k] = ref_ft.pos[i - 1];
      } else if (q[i] == 2) { /* We are at the end of a signal-rich region
      				 we are closing it, without creating a new one 
				(we're going into a signal-poor region) */
        Reg[0][k] = ref_ft.pos[i];
	cnt[k] += ref_ft.cnt[i];
	k++;
        if ((unsigned int)k >= mLen3) {
          mLen3 *= 2;
          for(j = 0; j < 2; j++) {
            if ((Reg[j] = (int*)realloc(Reg[j], mLen3 * sizeof(int))) == NULL) {
              perror("split_seq: realloc");
              exit(1);
	    }
          } 
        }
        Reg[0][k] = 0;
        Reg[1][k] = 0;
      } else { /* Extend signal-rich region */
        cnt[k] += ref_ft.cnt[i];
      }
    }
#ifdef DEBUG
    fprintf (stderr, "%d, %d, %d, %d, %d\n", i, k, Reg[0][k], Reg[1][k], cnt[k]);
#endif
  }
  /* Print out signal enriched regions */
  if (Reg[1][k] == 0) k--;
  if (Reg[0][k] == 0) Reg[0][k] = 1;
#ifdef DEBUG
  fprintf(stderr, "Printout number of Regions: k= %d\n", k);
#endif
  for (i = k; i >= 1; i--) {
    /* Print out (SGA format) */
    /* printf("%s\t%s\t%d\t%c\t%d\n", ref_ft.seq_id, ref_ft.name[i], Reg[0][i], '+', cnt[i]);
    printf("%s\t%s\t%d\t%c\t%d\n", ref_ft.seq_id, ref_ft.name[i], Reg[1][i], '-', cnt[i]);
    */
    printf("%s\t%s\t%d\t%c\t%d\n", ref_ft.seq_id, ref_ft.name[i], Reg[0][i], '+', cnt[i]);
    printf("%s\t%s\t%d\t%c\t%d\n", ref_ft.seq_id, ref_ft.name[i], Reg[1][i], '-', cnt[i]);
    seq_len += Reg[1][i] - Reg[0][i];
    tot_cnts += cnt[i];
  }
  if (end != 0)
    printf("%s\t%s\t%d\t%c\t%d\n", ref_ft.seq_id, ref_ft.name[len + 1], ref_ft.pos[len + 1], '0', 1);
  if (k && (len != 0)) {
    /*
    fprintf(stderr,"# %s :  Number of segments %d, Tot. length %d, Ave. length %d, Tot. counts %d\n", ref_ft.seq_id, k, seq_len, seq_len/k, tot_cnts); 
    */
    NumSegs += k;
    TotFragLen += seq_len;
    AvFragLen += seq_len/k;
    TotCntsFrag += tot_cnts;
    AvCntsFrag += tot_cnts/k;
  }
  if (last_seq) {
    fprintf(stderr,"# Num of Sequences : %d , Total Sequence Length : %lu (bp) , Total Counts : %lu , Total Num of Fragments : %d\n\n", TotSeqs, TotLen, TotCnts, NumSegs);
    fprintf(stderr,"# Total Fragment Length : %lu , Average Fragment Length : %lu (bp) , Percentage of Total Length : %2.5f\n\n", TotFragLen, AvFragLen/TotSeqs, (float)TotFragLen/(float)TotLen);
    fprintf(stderr,"# Percentage of Total Counts : %2.5f , Average Num of Counts per Fragment : %lu , Num of Counts per bp : %2.5f \n", (float)TotCntsFrag/(float)TotCnts, AvCntsFrag/TotSeqs, (float)TotCntsFrag/(float)TotFragLen);
  }
 
  free(q);
  free(cnt);
  for(j = 0; j < 2; j++) {
    if (Reg[j] != NULL)
      free(Reg[j]);
  } 
  free(Reg);
}

int
process_sga(FILE *input, char *iFile) 
{
  char seq_id_prev[SEQ_ID] = "";
  int pos, cnt, end = 0, last_pos = 0;
  size_t mLen = BUF_SIZE;
  char *s, *res, *buf;
  size_t bLen = LINE_SIZE;
  unsigned int k = 0;
  unsigned int i = 0;
  int end_flag = 0;

  if (options.debug && input != stdin) {
    char sort_cmd[1024] = "sort -s -c -k1,1 -k3,3n ";
    fprintf(stderr, "Check whether file %s is properly sorted\n", iFile);
    if (strcat(sort_cmd, iFile) == NULL) {
      fprintf(stderr, "strcat failed\n");
      return 1;
    }
    if (strcat(sort_cmd, " 2>/tmp/sortcheck.out") == NULL) {
      fprintf(stderr, "strcat failed\n");
      return 1;
    }
    fprintf(stderr, "executing : %s\n", sort_cmd);
    int sys_code = system(sort_cmd);
    /*fprintf(stderr,"system command sort : return code %d\n", sys_code);*/
    if (sys_code != 0) {
      fprintf(stderr, "system command failed\n");
      return 1;
    }
    struct stat file_status;
    if(stat("/tmp/sortcheck.out", &file_status) != 0){
      fprintf(stderr, "could not stat\n");
      return 1;
    }
    if (file_status.st_size != 0) {
      fprintf(stderr, "SGA file %s is not properly sorted\n", iFile);
      return 1;
    } else {
      system("/bin/rm /tmp/sortcheck.out");
    }
  }
  if ((ref_ft.pos = (int*)calloc(mLen, sizeof(int))) == NULL) {
    perror("process_sga: malloc");
    exit(1);
  }
  if ((ref_ft.strand = (char*)calloc(mLen, sizeof(char))) == NULL) {
    perror("process_sga: malloc");
    exit(1);
  }
  if ((ref_ft.cnt = (int*)calloc(mLen, sizeof(int))) == NULL) {
    perror("process_sga: malloc");
    exit(1);
  }
  if ((ref_ft.name = (char**)calloc(mLen, sizeof(*(ref_ft.name)))) == NULL) {
    perror("process_sga: malloc");
    exit(1);
  }
  if ((s = malloc(bLen * sizeof(char))) == NULL) {
    perror("process_sga: malloc");
    exit(1);
  }
#ifdef DEBUG
  int c = 1; 
#endif
/*
  while (fscanf(f,"%s %s %d %c %d", seq_id, ft, &pos, &strand, &cnt) != EOF) {
  */
  while ((res = fgets(s, (int) bLen, input)) != NULL) {
    size_t cLen = strlen(s);
    char seq_id[SEQ_ID] = "";
    char ft[FT_MAX] = ""; 
    char position[POS_MAX] = ""; 
    char count[CNT_MAX] = ""; 
    char strand = '\0';
    char ext[EXT_MAX]; 
    unsigned int j = 0;

    memset(ext, 0, (size_t)EXT_MAX);
    while (cLen + 1 == bLen && s[cLen - 1] != '\n') {
      bLen *= 2;
      if ((s = realloc(s, bLen)) == NULL) {
        perror("process_file: realloc");
        exit(1);
      }
      res = fgets(s + cLen, (int) (bLen - cLen), input);
      cLen = strlen(s);
    }
    if (s[cLen - 1] == '\n')
      s[cLen - 1] = 0;

    buf = s;
    /* Get SGA fields */
    /* SEQ ID */
    while (*buf != 0 && !isspace(*buf)) {
      if (j >= SEQ_ID) {
        fprintf(stderr, "Seq ID is too long \"%s\" \n", buf);
        exit(1);
      }
      seq_id[j++] = *buf++;
    }
    while (isspace(*buf))
      buf++;
    /* FEATURE */
    j = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (j >= FT_MAX) {
        fprintf(stderr, "Feature is too long \"%s\" \n", buf);
        exit(1);
      }
      ft[j++] = *buf++;
    }
    while (isspace(*buf))
      buf++;
    /* Position */
    j = 0;
    while (isdigit(*buf)) { 
      if (j >= POS_MAX) {
        fprintf(stderr, "Position is too large \"%s\" \n", buf);
        exit(1);
      }
      position[j++] = *buf++;
    }
    position[j] = 0;
    pos = atoi(position);
    while (isspace(*buf))
      buf++;
    /* Strand */
    strand = *buf++;
    while (isspace(*buf))
      buf++;
    /* Counts */
    j = 0;
    while (isdigit(*buf)) { 
      if (j >= CNT_MAX) {
        fprintf(stderr, "Count is too large \"%s\" \n", buf);
        exit(1);
      }
      count[j++] = *buf++;
    }
    count[j] = 0;
    cnt = atoi(count);
    while (isspace(*buf))
      buf++;
    /* SGA Extension */
    j = 0;
    while (*buf != 0) {
      if (j >= EXT_MAX) {
        fprintf(stderr, "Extension is too long \"%s\" \n", buf);
        exit(1);
      }
      ext[j++] = *buf++;
    }
#ifdef DEBUG
    printf(" [%d] seq ID: %s   Feat: %s (%c)  Pos: %d  Cnts: %d Ext: %s\n", c++, seq_id, ft, strand, pos, cnt, ext);
#endif
    if (k >= mLen - 1) {
      mLen *= 2;
#ifdef DEBUG
      fprintf(stderr, "reallocating memory for ref_ft.pos ref_ft.strand ref_ft.cnt (k=%d, size=%d)\n", c, mLen);
#endif
      if ((ref_ft.pos = (int *)realloc(ref_ft.pos, mLen * sizeof(int))) == NULL) {
    	perror("process_sga: realloc");
    	exit(1);
      }
      if ((ref_ft.cnt = (int *)realloc(ref_ft.cnt, mLen * sizeof(int))) == NULL) {
    	perror("process_sga: realloc");
    	exit(1);
      }
      if ((ref_ft.strand = (char *)realloc(ref_ft.strand, mLen * sizeof(char))) == NULL) {
    	perror("process_sga: realloc");
    	exit(1);
      }
      if ((ref_ft.name = (char**)realloc(ref_ft.name, mLen * sizeof(*(ref_ft.name)))) == NULL) {
        perror("process_sga: malloc");
        exit(1);
      }
    }
    /* Check Chromosome BEGINNING, split sequences in regions and printout results*/
    if (strcmp(seq_id, seq_id_prev) != 0) {
      split_seq((int)k, end, 0); 
      for (i = 1; i <= k; i++) { 
	if (ref_ft.name[i] != NULL)
          free(ref_ft.name[i]);
      }
      if (end_flag) {
	if (ref_ft.name[k + 1] != NULL)
          free(ref_ft.name[k + 1]);
      }
      k = 0;
      end = 0;
      end_flag = 0;
      TotLen += last_pos;
      TotSeqs ++;
      ref_ft.pos[0] = 1;
      ref_ft.cnt[0] = 0;
      strcpy(seq_id_prev, seq_id);
      strcpy(ref_ft.seq_id, seq_id);
    }
    if (ref_ft.ft == NULL) {
      k++;
      ref_ft.name[k] = malloc(strlen(ft) + 1);
      strcpy(ref_ft.name[k], ft);
      ref_ft.strand[k] = strand;
      ref_ft.pos[k] = pos;
      if (cnt > Coff)
          ref_ft.cnt[k] = Coff;
      else
          ref_ft.cnt[k] = cnt;
      TotCnts += ref_ft.cnt[k];
    } else if (ref_ft.ft_str == '\0') {
      if (strcmp(ft, ref_ft.ft) == 0) {
        k++;
        ref_ft.name[k] = malloc(strlen(ft) + 1);
        strcpy(ref_ft.name[k], ft);
        ref_ft.strand[k] = strand;
        ref_ft.pos[k] = pos;
	if (cnt > Coff)
          ref_ft.cnt[k] = Coff;
        else
          ref_ft.cnt[k] = cnt;
      TotCnts += ref_ft.cnt[k];
      }
    } else if (strand_flag == 1) {
      if (strand == ref_ft.ft_str) {
        k++;
        ref_ft.name[k] = malloc(strlen(ft) + 1);
        strcpy(ref_ft.name[k], ft);
        ref_ft.strand[k] = strand;
        ref_ft.pos[k] = pos;
	if (cnt > Coff)
          ref_ft.cnt[k] = Coff;
        else
          ref_ft.cnt[k] = cnt;
      TotCnts += ref_ft.cnt[k];
      }
    } else if (strcmp(ft, ref_ft.ft) == 0  && strand == ref_ft.ft_str) {
        k++;
        ref_ft.name[k] = malloc(strlen(ft) + 1);
        strcpy(ref_ft.name[k], ft);
        ref_ft.strand[k] = strand;
        ref_ft.pos[k] = pos;
	if (cnt > Coff)
          ref_ft.cnt[k] = Coff;
        else
          ref_ft.cnt[k] = cnt;
      TotCnts += ref_ft.cnt[k];
    } 
    if (strcmp(ft, "END") == 0 ) {
        ref_ft.name[k + 1] = malloc(strlen(ft) + 1);
        strcpy(ref_ft.name[k + 1], ft);
        end = pos;
        end_flag = 1;
    }
    last_pos = pos;
#ifdef DEBUG
    fprintf(stderr,"k = [%d]   pos = %d cnts = %d\n", k, ref_ft.pos[k], ref_ft.cnt[k]);
#endif
  } /* End of While */
  free(s);
  /* Partition last chromosome */ 
  TotLen += last_pos;
  split_seq((int)k, end, 1); 
  for (i = 1; i <= k; i++) { 
    if (ref_ft.name[i] != NULL)
      free(ref_ft.name[i]);
  }
  if (end_flag) {
    if (ref_ft.name[k + 1] != NULL)
      free(ref_ft.name[k + 1]);
  }
  free(ref_ft.name);
  free(ref_ft.pos);
  free(ref_ft.cnt);
  free(ref_ft.strand);
  if (input != stdin) {
    fclose(input);
  }
  return 0;
}

int
main(int argc, char *argv[])
{
#ifdef DEBUG
  mcheck(NULL);
  mtrace();
#endif
  FILE *input;

  while (1) {
    int c = getopt(argc, argv, "f:dhs:p:c:");
    if (c == -1)
      break;
    switch (c) {
      case 'f':
        Feature = optarg;
        break;
      case 'd':
        options.debug = 1;
        break;
      case 'h':
        options.help = 1;
        break;
    case 's':
        Dthres = atof(optarg);
	break;
    case 'p':
        Tpen = atoi(optarg);
	break;
    case 'c':
        Coff = atoi(optarg);
	break;
      default:
        printf ("?? getopt returned character code 0%o ??\n", c);
    }
  }
  if (optind > argc || options.help == 1 || Dthres == 0 || Tpen == 0 || Coff < 0) {
    fprintf(stderr, "Usage: %s [options] [-f <feature>] -s <densitiy threshold> -p <transition penalty> [<] <SGA file>\n"
             "      - version %s\n"
             "      where options are:\n"
	     "  \t\t -h     Show this help text\n"
	     "  \t\t -d     Print debug information and check SGA file\n"
	     "  \t\t -c     Count Cut-off (default is %d). It must be >= 0.\n"
	     "\n\tPartitioning Tool.\n"
	     "\n\tThe program reads a ChIP-seq data file (or from stdin [<]) in SGA format (<SGA file>),\n"
	     "\tand heuristically partitions the data corresponding to a specific feature\n"
	     "\t<feature>, if the latter is given.\n"
	     "\tThe <feature> parameter is a name that corresponds to the second field of the SGA file.\n"
	     "\tIt might optionally include the strand specification (+|-).\n"
	     "\tIf no feature is given then all input tags are processed.\n"
	     "\tThe SGA input file MUST BE sorted by sequence name (or chromosome id), position,\n"
	     "\tand strand.\n"
	     "\tOne should check the input SGA file with the following command:\n"
	     "\tsort -s -c -k1,1 -k3,3n -k4,4 <SGA file>.\n\n"
	     "\tIn debug mode (-d), the program performs the sorting order check.\n\n"
	     "\tInput parameters are the density threshold (<density threshold>),\n"
	     "\tthe transition penalty (<transition penalty>), which are used to\n"
             "\tsplit the data in regions of high density count.\n" 
	     "\tA value can be optionally specified as a cut-off for the feature counts.\n"
             "\tThe output is an SGA-formatted list containing the regions of interest.\n"
             "\tThe program also generates (to stderr) a statistical report with the\n"
             "\tfollowing information:\n"
	     "\t  - Total number of processed sequences, total DNA length, and\n" 
	     "\t    total number of fragments;\n"
	     "\t  - Total length of fragments, average fragment length, and\n"
	     "\t    percentage of total DNA length;\n"
             "\t  - Percentage of total counts, average number of counts, and\n"
	     "\t    number of count per bp (count density).\n\n",
	     argv[0], VERSION, Coff);
      return 1;
  }
  if (argc > optind) {
      if(!strcmp(argv[optind],"-")) {
          input = stdin;
      } else {
          input = fopen(argv[optind], "r");
          if (NULL == input) {
              fprintf(stderr, "Unable to open '%s': %s(%d)\n",
                  argv[optind], strerror(errno), errno);
             exit(EXIT_FAILURE);
          }
          if (options.debug)
             fprintf(stderr, "Processing file %s\n", argv[optind]);
      }
  } else {
      input = stdin;
  }
  if (options.debug) {
    fprintf(stderr, " Arguments:\n");
    fprintf(stderr, " Selected Feature : %s\n", Feature);
    fprintf(stderr, " Density Threshold : %f\n\n", Dthres);
    fprintf(stderr, " Transition Penalty : %d\n\n", Tpen);
    fprintf(stderr, " Count Cut-off : %d\n\n", Coff);
  }
  /* Process Feature Specs */
  if (Feature == NULL) {
    ref_ft.ft = NULL;      /* Process all features */
    ref_ft.ft_str = '\0';
  } else {
    ref_ft.ft = malloc(FT_MAX * sizeof(char));
    char *s = Feature;
    int i = 0;
    while (*s != 0  && !isspace(*s)) {
      if (i >= FT_MAX) {
        fprintf(stderr, "Feature Description too long \"%s\" \n", Feature);
        return 1;
      }
      ref_ft.ft[i++] = *s++;
    }
    ref_ft.ft[i] = '\0';
    ref_ft.ft_str = '\0';
    while (isspace(*s++))
    ref_ft.ft_str = *s;
  }
  if (options.debug) {
    if (ref_ft.ft_str == '\0' && ref_ft.ft == NULL) {
      fprintf(stderr, "Feature Specs: ALL -> Process all features\n");
    } else if (ref_ft.ft_str == '\0') {
      fprintf(stderr, "Feature Specs: Feature name : %s\n", ref_ft.ft);
    } else {
      fprintf(stderr, "Feature Specs: Feature name/str : %s %c\n", ref_ft.ft, ref_ft.ft_str);
    }
  }
  if ( ref_ft.ft != NULL && (strcmp(ref_ft.ft, "+") == 0 || strcmp(ref_ft.ft, "-") == 0)) {
    strcpy(&ref_ft.ft_str, ref_ft.ft);
    strand_flag = 1;
    if (options.debug)
      fprintf(stderr, "Feature Specs: Process all features on str : %c\n", ref_ft.ft_str);
  }
  if (process_sga(input, argv[optind++]) != 0) {
    return 1;
  }
  free(ref_ft.ft);
  return 0;
}