/*
  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 SIB 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"
#include "sysck.h"

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

/* Some versions of GCC won't warn about stdlib's system()'s result being
 * ignored unless we give it this attribute. This is the case even if we have
 * all the -Wall, -Wextra etc flags switched on. This allows me to reproduce the
 * bug mentioned by SM on 2023-12-01 (see Teams chat).*/
// TODO: this will become moot anyway as soon as we stop using tmp files by
// checking the exit status of system(3) rather than the STDERR of sort(1).
int system(const char *command) __attribute__((warn_unused_result));

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 -k4,4 ";
        /*fprintf(stderr, "Check whether file %s is properly sorted\n",
                iFile);*/
        if (strcat(sort_cmd, iFile) == NULL) {
            fprintf(stderr, "strcat failed\n");
            return 1;
        }
        int sys_code = system(sort_cmd);
        if (sys_code != 0) {
            fprintf(stderr,"system command %s : return code %d\n", sort_cmd, sys_code);
            fprintf(stderr, "sorting command command failed\n");
            return 1;
        }
    }
    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:dht:s: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':
            fprintf(stderr, "WARNING: option -s is deprecated; please use option -t instead.\n");
            Dthres = atof(optarg);
            break;
        case 't':
            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>] -t <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 -s          Deprecated, please use -t\n"
                "  \t\t -d          Print debug information and check SGA file\n"
                "  \t\t -c <cutoff> 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;
}
