/*
  chip_peak.c

  Signal Peak Tool.
  The program locates signal peaks within a SGA file.

  # Arguments:
  # feature type, strand, integration range (Win1)
  # minimal distance (Win2), counts threshold (Thres)

  Giovanna Ambrosini, EPFL/ISREC, 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>
#include "hashtable.h"
#ifdef DEBUG
#include <mcheck.h>
#endif

#include "version.h"

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

typedef struct _options_t {
    int help;
    int debug;
    int db;
    char *dbPath;
    int refine;
    int oriented;
} options_t;

static options_t options;

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

feature_t ref_ft;
feature_t ref_ft_plus;
feature_t ref_ft_minus;

typedef struct _locmax_t {
    unsigned long *pos;
    int *cnt;
    char **name;
} locmax_t, *locmax_p_t;

locmax_t lm_plus;
locmax_t lm_minus;

int strand_flag = 0;

static hash_table_t *size_table = NULL;

char *Feature = NULL;
int Win1 = 0;
int Win2 = 0;
int Thres = 50;
int Coff = 1;
char *progname = NULL;

void show_help(void) {
    fprintf(stderr,
            "Usage: %s [options] [-f <feature>] -t <threshold> -w <window> -v <vicinity> [<] <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 -i <dir>   Directory where the chr_size file is found\n"
            "  \t\t            [Default=$CWD]\n"
            "  \t\t -o         Oriented strand processing\n"
            "  \t\t -r         Refine Peak Positions\n"
            "  \t\t -c         Count Cut-off (default is %d)\n"
            "\n\tLocate signal peaks within SGA files.\n"
            "\n\tThe program reads a ChIP-seq data file (or from stdin [<]) in SGA format (<SGA file>),\n"
            "\tand detects signal peaks for ChIP-tag positions corresponding to a specific genomic\n"
            "\tfeature (<feature>). The <feature> parameter is a name that corresponds to the second \n"
            "\tfield of the input SGA file. It might optionally include the strand specification (+|-).\n"
            "\tIf no feature name is given then all input tags are processed.\n"
            "\tIf the oriented option is specified (-o), peaks are separately detected on plus and minus\n"
            "\tstrands, respectively.\n"
            "\tThe SGA input file MUST BE sorted by sequence name (or chromosome id), position, and 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"
            "\tAdditional input parameters are the integration range (<window>), the minimal distance\n"
            "\tamongst a group of high count local peaks (<vicinity>), and the peak threshold (<threshold>)\n"
            "\twhose default is %d. A value can be optionally specified as a cut-off for the feature counts.\n"
            "\tIf peak refinement is required (-r option), the program recomputes the position of each\n"
            "\tpeak by taking the center of gravity of the read counts within the peak region (<window>).\n"
            "\tThe program also checks whether, after peak refinement, the new peak coordinates are still\n"
            "\twithin chromosome boundaries. For that reason, the file chr_size must be read.\n"
            "\tThe output is an SGA-formatted list containing signal peak locations.\n\n",
            progname, VERSION, Coff, Thres);
}

int process_size(void) {
    FILE *input;
    //int c;
    char buf[LINE_SIZE];
    char *chrSizeFile;
    int cLen;

    if (options.db) {
        cLen = (int) strlen(options.dbPath) + 10;
        if ((chrSizeFile = (char *) malloc(cLen * sizeof(char))) == NULL) {
            perror("process_ac: malloc");
            exit(1);
        }
        strcpy(chrSizeFile, options.dbPath);
    } else {
        cLen = 16 + 10;
        if ((chrSizeFile = (char *) malloc(cLen * sizeof(char))) == NULL) {
            perror("process_ac: malloc");
            exit(1);
        }
        /*
         * Replace /local/db/genome (which cannot be assumed to exist on users'
         * machines) with $CWD.
         strcpy(chrSizeFile, "/local/db/genome");
         */
        strcpy(chrSizeFile, ".");
    }
    strcat(chrSizeFile, "/chr_size");

    input = fopen(chrSizeFile, "r");
    if (input == NULL) {
        fprintf(stderr,
                "ERROR: Could not open file %s - %s. See option -i.\n\n",
                chrSizeFile, strerror(errno));
        show_help();
        return 1;
    }

    size_table = hash_table_new(MODE_COPY);
    while (fgets(buf, LINE_SIZE, input) != NULL) {
        char *s;
        char chr_size[CHR_SIZE] = "";
        char ncbi_ac[AC_MAX] = "";
        int i = 0;
        int ac_len = 0;
        int size_len = 0;
        /*int valid = 1; */
        s = buf;
        /* Check line */
        /* Get first character: if # skip line */
        if (*s == '#')
            continue;
        /* Chromosome NCBI AC */
        i = 0;
        while (*s != 0 && !isspace(*s)) {
            if (i >= AC_MAX) {
                fprintf(stderr, "AC too long \"%s\" \n", s);
                fclose(input);
                exit(1);
            }
            ncbi_ac[i++] = *s++;
        }
        if (i < AC_MAX)
            ncbi_ac[i] = 0;
        ac_len = i + 1;
        while (isspace(*s))
            s++;
        i = 0;
        /* Chrom SIZE */
        while (*s != 0 && !isspace(*s)) {
            if (i >= CHR_SIZE) {
                fprintf(stderr, "Size too long in %s\n", s);
                fclose(input);
                exit(1);
            }
            chr_size[i++] = *s++;
        }
        if (i < CHR_NB)
            chr_size[i] = 0;
        size_len = i + 1;

        /* printf("adding key %s (len = %d) value %s (ac) (len = %d) to hash table\n", ncbi_ac, ac_len, chr_size, size_len); */
        hash_table_add(size_table, ncbi_ac, (size_t) ac_len, chr_size,
                       (size_t) size_len);
        if (options.debug) {
            char *csize =
                hash_table_lookup(size_table, ncbi_ac, (size_t) ac_len);
            fprintf(stderr,
                    " SIZE Hash table: %s (len = %d) -> %s (len = %d)\n",
                    ncbi_ac, ac_len, csize, size_len);
        }
    }
    return 0;
}

void locate_peaks(int len) {
    /* Compute sum in window Win1; store high values in shorter arrays (npo, nct) */
    unsigned long long sum = 0;
    int cnts = 0;
    int i, k;
    int j = 0;
    size_t mLen = BUF_SIZE;
    unsigned int size;
    int *lm;
    char ft_name[FT_MAX + 4];
    char str = '0';
    int ac_len = 0;
    char *csize = NULL;
    int chr_size = 0;

    if (ref_ft.ft_str != '\0') {
        str = ref_ft.ft_str;
    }
    if (ref_ft.ft != NULL) {
        strcpy(ft_name, ref_ft.ft);
        strcat(ft_name, "_p");
    }
    if (options.refine) {
        /* Get Chromosome size */
        ac_len = (int) strlen(ref_ft.seq_id) + 1;
        csize =
            hash_table_lookup(size_table, ref_ft.seq_id, (size_t) ac_len);
        if (csize != NULL) {
            chr_size = (int) atoi(csize);
        }
    }
    for (i = 1; i <= len; i++) {
        sum = ref_ft.cnt[i];
        /* Sum up all tag counts within the window range +-Win1/2 */
        for (k = i - 1;
             (long) ref_ft.pos[k] >= (long) (ref_ft.pos[i] - Win1 / 2);
             k--) {
            sum += ref_ft.cnt[k];
        }
        for (k = i + 1; ref_ft.pos[k] <= ref_ft.pos[i] + Win1 / 2; k++) {
            sum += ref_ft.cnt[k];
        }
        if ((unsigned int) j >= mLen - 1) {
            mLen *= 2;
            if ((ref_ft.npo =
                 (unsigned long *) realloc(ref_ft.npo,
                                           mLen *
                                           sizeof(unsigned long))) ==
                NULL) {
                perror("locate_peaks: realloc");
                exit(1);
            }
            if ((ref_ft.nct =
                 (int *) realloc(ref_ft.nct,
                                 mLen * sizeof(int))) == NULL) {
                perror("locate_peaks: realloc");
                exit(1);
            }
            if ((ref_ft.ptr =
                 (int *) realloc(ref_ft.ptr,
                                 mLen * sizeof(int))) == NULL) {
                perror("locate_peaks: realloc");
                exit(1);
            }
        }
        if (sum >= (unsigned long long) Thres) {
            j++;
            ref_ft.npo[j] = ref_ft.pos[i];
            ref_ft.nct[j] = sum;
            ref_ft.ptr[j] = i;
        }
    }
    /* Initialize Local Maxima Array lm */
    size = (unsigned int) j + 1;
    if ((lm = (int *) calloc((size_t) size, sizeof(size_t))) == NULL) {
        perror("locate_peaks: calloc");
        exit(1);
    }
    /* Keep only one maximum value (peak) within a vicinity range +-Win2 */
    /* Record local maxima in lm flag Array */
    /* We distinguish three different cases :
       - local maxima within Win2 distance in forward direction (lm = 1)
       - local maxima within Win2 distance in backward direction (lm = 2)
       - local maxima within +-Win2 in both forw/back directions (lm = 3)
     */
    /* Select maxima forward path (alike segmentation algorithm) */
    int max = 1;
    for (i = 2; i <= j; i++) {
        if (ref_ft.npo[i] > ref_ft.npo[max] + Win2) {
            /* if the distance between two local
               maxima (i, max) is greater than Win2
               keep pos max as a local maxima and
               increment lm flag */
            lm[max]++;
            max = i;
        } else if (ref_ft.nct[i] > ref_ft.nct[max]) {   /* Else, max is not a local maxima */
            max = i;
        }
    }
    lm[max]++;
    /* Select maxima backward path */
    max = j;
    for (i = j - 1; i >= 0; i--) {
        if (ref_ft.npo[i] < ref_ft.npo[max] - Win2) {
            /* if the distance between two local
               maxima (i, max) is greater than Win2
               keep pos max as a local maxima and
               increment by 2 lm flag */
            lm[max] += 2;
            max = i;
        } else if (ref_ft.nct[i] >= ref_ft.nct[max]) {  /* Else, max is not a local maxima */
            max = i;
        }
    }
    lm[max] += 2;

    /* Print out local maxima positions */
    /* Only positions with lm[i]=3 should be considered as signal peaks */
    for (i = 1; i <= j; i++) {
#ifdef DEBUG
        if (ref_ft.ft != NULL) {
            fprintf(stderr, "dbg: %s\t%s\t%lu\t%c\t%d\t%d\n",
                    ref_ft.seq_id, ft_name, ref_ft.npo[i], str,
                    ref_ft.nct[i], lm[i]);
        } else {
            fprintf(stderr, "dbg: %s\t%s\t%lu\t%c\t%d\t%d\n",
                    ref_ft.seq_id, ref_ft.name[i], ref_ft.npo[i], str,
                    ref_ft.nct[i], lm[i]);
        }
#endif
        if (lm[i] == 3) {
            if (!options.refine) {
                if (ref_ft.ft != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", ref_ft.seq_id, ft_name,
                           ref_ft.npo[i], str, ref_ft.nct[i]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", ref_ft.seq_id,
                           ref_ft.name[i], ref_ft.npo[i], str,
                           ref_ft.nct[i]);
                }
            } else {
                double sum2;
                /* Refine peak positions */
                /*
                   printf("\ndbg: Refine peak pos for POS %d i=%d : \tk=%d : ref_ft.cnt %d\tref_ft.pos  %d\n", ref_ft.npo[i], i, ref_ft.ptr[i], ref_ft.cnt[ref_ft.ptr[i]], ref_ft.pos[ref_ft.ptr[i]]);
                 */
                sum2 =
                    (double) (ref_ft.npo[i] * ref_ft.cnt[ref_ft.ptr[i]]);
                cnts = ref_ft.cnt[ref_ft.ptr[i]];
                /*
                   printf("\ndbg:  INIT: sum %llu   cnts %d \n\n", sum, cnts);
                 */
                /* Recompute peak position within the window range +-(Win1/2) */
                for (k = ref_ft.ptr[i] - 1;
                     ref_ft.pos[k] >= ref_ft.npo[i] - Win1 / 2; k--) {
                    /*
                       printf("i=%d ref_ft.npo %d\tk=%d : ref_ft.cnt %d\tref_ft.pos  %d\n", i, ref_ft.npo[i], k, ref_ft.cnt[k], ref_ft.pos[k]);
                     */
                    sum2 += (double) (ref_ft.cnt[k] * ref_ft.pos[k]);
                    cnts += ref_ft.cnt[k];
                }
                for (k = ref_ft.ptr[i] + 1;
                     ref_ft.pos[k] <= ref_ft.npo[i] + Win1 / 2; k++) {
                    /*
                       printf("i=%d ref_ft.npo %d\tk=%d : ref_ft.cnt %d\tref_ft.pos  %d\n", i, ref_ft.npo[i], k, ref_ft.cnt[k], ref_ft.pos[k]);
                     */
                    sum2 += (double) (ref_ft.cnt[k] * ref_ft.pos[k]);
                    cnts += ref_ft.cnt[k];
                }
                /*
                   printf("\ndbg: OLD POS %d : \t", ref_ft.npo[i]);
                   printf("SUM = %llu CNTS = %d NEW POS : %int \n\n", sum2, cnts, (int)(sum2/cnts));
                 */
                ref_ft.npo[i] = (unsigned long) (sum2 / cnts);
                /* Check Chromosome Boundaries */
                if (chr_size) {
                    if (ref_ft.npo[i] > (unsigned int) chr_size) {
                        fprintf(stderr,
                                "WARNING: peak position %s\t%lut%c\t%d goes beyond chromosome boundaries (chrom size = %d)\n",
                                ref_ft.seq_id, ref_ft.npo[i], str,
                                ref_ft.nct[i], chr_size);
                        continue;
                    }
                }
                if (ref_ft.ft != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", ref_ft.seq_id, ft_name,
                           ref_ft.npo[i], str, ref_ft.nct[i]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", ref_ft.seq_id,
                           ref_ft.name[i], ref_ft.npo[i], str,
                           ref_ft.nct[i]);
                }
            }                   /* If refine peak pos */
        }
    }
    free(lm);
}

int
locate_peaks_oriented(feature_p_t feat, locmax_p_t lmax, char strand,
                      int len) {
    /* Compute sum in window Win1; store high values in shorter arrays (npo, nct) */
    unsigned long long sum = 0;
    int cnts = 0;
    int i, k = 0;
    int j = 0;
    size_t mLen = BUF_SIZE;
    unsigned int size;
    int *lm;
    char ft_name[FT_MAX + 4];
    int ac_len = 0;
    char *csize = NULL;
    int chr_size = 0;

    /* Feature name and strand are only passed for debugging reasons: they can be omitted */
    if (ref_ft.ft != NULL) {
        strcpy(ft_name, ref_ft.ft);
        strcat(ft_name, "_p");
    }
    if (options.refine) {
        /* Get Chromosome size */
        ac_len = (int) strlen(feat->seq_id) + 1;
        csize =
            hash_table_lookup(size_table, feat->seq_id, (size_t) ac_len);
        if (csize != NULL) {
            chr_size = (int) atoi(csize);
        }
    }
    /*printf("locate_peaks_oriented : strand %c len %d feature %s Win %d\n", strand, len, ft_name, Win1); */
    /*
       for (i = 1; i <=10; i++) {
       printf ("%s ft pos: %lu  ft cnts: %d\n", ft_name, feat->pos[i], feat->cnt[i]);
       } */
    for (i = 1; i <= len; i++) {        /* Loop on Feature Array */
        sum = feat->cnt[i];
        /* Sum up all tag counts within the window range +-Win1/2 */
        for (k = i - 1;
             (long) feat->pos[k] >= (long) (feat->pos[i] - Win1 / 2);
             k--) {
            sum += feat->cnt[k];
        }
        for (k = i + 1; feat->pos[k] <= feat->pos[i] + Win1 / 2; k++) {
            sum += feat->cnt[k];
        }
        if ((unsigned int) j >= mLen - 1) {
            mLen *= 2;
            if ((feat->npo =
                 (unsigned long *) realloc(feat->npo,
                                           mLen *
                                           sizeof(unsigned long))) ==
                NULL) {
                perror("locate_peaks_oriented: realloc");
                exit(1);
            }
            if ((feat->nct =
                 (int *) realloc(feat->nct, mLen * sizeof(int))) == NULL) {
                perror("locate_peaks_oriented: realloc");
                exit(1);
            }
            if ((feat->ptr =
                 (int *) realloc(feat->ptr, mLen * sizeof(int))) == NULL) {
                perror("locate_peaks_oriented: realloc");
                exit(1);
            }
        }
        /*printf ("feature %d : pos: %lu j: %d SUM: %llu Thres: %d\n", i, feat->pos[i], j, sum, Thres); */
        if (sum >= (unsigned long long) Thres) {
            j++;
            feat->npo[j] = feat->pos[i];
            feat->nct[j] = sum;
            feat->ptr[j] = i;
        }
    }
    /* Initialize Local Maxima Array lm */
    size = (unsigned int) j + 1;
    if ((lm = (int *) calloc((size_t) size, sizeof(size_t))) == NULL) {
        perror("locate_peaks: calloc");
        exit(1);
    }
    /* Keep only one maximum value (peak) within a vicinity range +-Win2 */
    /* Record local maxima in lm flag Array */
    /* We distinguish three different cases :
       - local maxima within Win2 distance in forward direction (lm = 1)
       - local maxima within Win2 distance in backward direction (lm = 2)
       - local maxima within +-Win2 in both forw/back directions (lm = 3)
     */
    /* Select maxima forward path (alike segmentation algorithm) */
    int max = 1;
    for (i = 2; i <= j; i++) {
        if (feat->npo[i] > feat->npo[max] + Win2) {
            /* If the distance between two local
               maxima (i, max) is greater than Win2
               keep pos max as a local maxima and
               increment lm flag */
            lm[max]++;
            max = i;
        } else if (feat->nct[i] > feat->nct[max]) {     /* Else, max is not a local maxima */
            max = i;
        }
    }
    lm[max]++;
    /* Select maxima backward path */
    max = j;
    for (i = j - 1; i >= 0; i--) {
        if (feat->npo[i] < feat->npo[max] - Win2) {
            /* If the distance between two local
               maxima (i, max) is greater than Win2
               keep pos max as a local maxima and
               increment by 2 lm flag */
            lm[max] += 2;
            max = i;
        } else if (feat->nct[i] >= feat->nct[max]) {    /* Else, max is not a local maxima */
            max = i;
        }
    }
    lm[max] += 2;

    /* Store local maxima positions into locmax Array */
    /* Only positions with lm[i]=3 should be considered as signal peaks */
    unsigned int n = 0;
    mLen = BUF_SIZE;
    for (i = 1; i <= j; i++) {
#ifdef DEBUG2
        if (ref_ft.ft != NULL) {
            fprintf(stderr, "dbg: %s\t%s\t%lu\t%c\t%d\t%d\n", feat->seq_id,
                    ft_name, feat->npo[i], strand, feat->nct[i], lm[i]);
        } else {
            fprintf(stderr, "dbg: %s\t%s\t%lu\t%c\t%d\t%d\n", feat->seq_id,
                    feat->name[i], feat->npo[i], strand, feat->nct[i],
                    lm[i]);
        }
#endif
        if (n >= mLen - 1) {
            mLen *= 2;
            if ((lmax->pos =
                 (unsigned long *) realloc(lmax->pos,
                                           mLen *
                                           sizeof(unsigned long))) ==
                NULL) {
                perror("locate_peaks_oriented: realloc");
                exit(1);
            }
            if ((lmax->cnt =
                 (int *) realloc(lmax->cnt, mLen * sizeof(int))) == NULL) {
                perror("locate_peaks_oriented: realloc");
                exit(1);
            }
            if (ref_ft.ft == NULL) {
                if ((lmax->name =
                     (char **) realloc(lmax->name,
                                       mLen * sizeof(*(lmax->name)))) ==
                    NULL) {
                    perror("process_sga: malloc");
                    exit(1);
                }
            }
        }
        if (lm[i] == 3) {
            if (!options.refine) {
                if (options.debug) {
                    if (ref_ft.ft != NULL) {
                        printf("%s\t%lu\t%c\t%d\n", ft_name, feat->npo[i],
                               strand, feat->nct[i]);
                    } else {
                        printf("%s\t%lu\t%c\t%d\n", feat->name[i],
                               feat->npo[i], strand, feat->nct[i]);
                    }
                }
                lmax->pos[n] = feat->npo[i];
                lmax->cnt[n] = feat->nct[i];
                if (ref_ft.ft == NULL) {
                    lmax->name[n] = malloc(strlen(feat->name[i]) + 1);
                    strcpy(lmax->name[n], feat->name[i]);
                }
                n++;
            } else {
                double sum2;
                /* Refine peak positions */
                sum2 = (double) (feat->npo[i] * feat->cnt[feat->ptr[i]]);
                cnts = feat->cnt[feat->ptr[i]];
                /* printf("\ndbg:  INIT: sum %llu   cnts %d \n\n", sum, cnts); */
                /* Recompute peak position within the window range +-(Win1/2) */
                for (k = feat->ptr[i] - 1;
                     feat->pos[k] >= feat->npo[i] - Win1 / 2; k--) {
                    sum2 += (double) (feat->cnt[k] * feat->pos[k]);
                    cnts += feat->cnt[k];
                }
                for (k = feat->ptr[i] + 1;
                     feat->pos[k] <= feat->npo[i] + Win1 / 2; k++) {
                    sum2 += (double) (feat->cnt[k] * feat->pos[k]);
                    cnts += feat->cnt[k];
                }
                /*
                   printf("\ndbg: OLD POS %d : \t", feat->npo[i]);
                   printf("SUM = %llu CNTS = %d NEW POS : %int \n\n", sum2, cnts, (int)(sum2/cnts));
                 */
                feat->npo[i] = (unsigned long) (sum2 / cnts);
                /* Check Chromosome Boundaries */
                if (chr_size) {
                    if (feat->npo[i] > (unsigned int) chr_size) {
                        fprintf(stderr,
                                "WARNING: peak position %s\t%lu\t%c\t%d goes beyond chromosome boundaries (chrom size = %d)\n",
                                feat->seq_id, feat->npo[i], strand,
                                feat->nct[i], chr_size);
                        continue;
                    }
                }
                if (options.debug) {
                    if (ref_ft.ft != NULL) {
                        printf("%s\t%lu\t%c\t%d\n", ft_name, feat->npo[i],
                               strand, feat->nct[i]);
                    } else {
                        printf("%s\t%lu\t%c\t%d\n", feat->name[i],
                               feat->npo[i], strand, feat->nct[i]);
                    }
                }
                lmax->pos[n] = feat->npo[i];
                lmax->cnt[n] = feat->nct[i];
                if (ref_ft.ft == NULL) {
                    lmax->name[n] = malloc(strlen(feat->name[i]) + 1);
                    strcpy(lmax->name[n], feat->name[i]);
                }
                n++;
            }                   /* If refine peak pos */
        }
    }
    free(lm);
    return n;
}

void merge_peaks(int m, int n, char *seq_id, char *ftname) {
    char ft_name[FT_MAX + 4];
    int i, k, j;

    if (ftname != NULL) {
        strcpy(ft_name, ftname);
        strcat(ft_name, "_p");
    }
    /*printf ("merge_peaks: m %d, n %d, seq_id %s , feature %s\n", m, n, seq_id, ft_name); */
    k = j = 0;
    for (i = 0; i < m + n;) {
        if (k < m && j < n) {
            if (lm_plus.pos[k] < lm_minus.pos[j]) {
                if (ftname != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id, ft_name,
                           lm_plus.pos[k], '+', lm_plus.cnt[k]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id,
                           lm_plus.name[k], lm_plus.pos[k], '+',
                           lm_plus.cnt[k]);
                }
                k++;
            } else {
                if (ftname != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id, ft_name,
                           lm_minus.pos[j], '-', lm_minus.cnt[j]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id,
                           lm_minus.name[j], lm_minus.pos[j], '-',
                           lm_minus.cnt[j]);
                }
                j++;
            }
            i++;
        } else if (k == m) {    /* print lm_minus array */
            for (; i < m + n;) {
                if (ftname != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id, ft_name,
                           lm_minus.pos[j], '-', lm_minus.cnt[j]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id,
                           lm_minus.name[j], lm_minus.pos[j], '-',
                           lm_minus.cnt[j]);
                }
                j++;
                i++;
            }
        } else {                /* j = n print lm_plus array */
            for (; i < m + n;) {
                if (ftname != NULL) {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id, ft_name,
                           lm_plus.pos[k], '+', lm_plus.cnt[k]);
                } else {
                    printf("%s\t%s\t%lu\t%c\t%d\n", seq_id,
                           lm_plus.name[k], lm_plus.pos[k], '+',
                           lm_plus.cnt[k]);
                }
                k++;
                i++;
            }
        }
    }
}

int process_sga(FILE *input, char *iFile) {
    char seq_id_prev[SEQ_ID] = "";
    int pos, cnt;
    size_t mLen = BUF_SIZE;
    size_t mLen1 = BUF_SIZE;
    char *s, *res, *buf;
    size_t bLen = LINE_SIZE;
    unsigned int k = 0;
    unsigned int j = 0;
    int k1 = 0, j1 = 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 (!options.oriented) {
        if ((ref_ft.pos =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft.cnt = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        /* Peak Array (position, count, feature name) */
        if ((ref_ft.npo =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft.nct = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft.ptr = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if (ref_ft.ft == NULL) {
            if ((ref_ft.name =
                 (char **) calloc(mLen, sizeof(*(ref_ft.name)))) == NULL) {
                perror("process_sga: malloc");
                exit(1);
            }
        }
    } else {
        /* Peak Arrays on both strands  (position, count, feature name) */
        if ((ref_ft_plus.pos =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_plus.cnt = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_plus.npo =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_plus.nct = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_plus.ptr = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if (ref_ft.ft == NULL) {
            if ((ref_ft_plus.name =
                 (char **) calloc(mLen,
                                  sizeof(*(ref_ft_plus.name)))) == NULL) {
                perror("process_sga: malloc");
                exit(1);
            }
        }
        if ((ref_ft_minus.pos =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_minus.cnt = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_minus.npo =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_minus.nct = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((ref_ft_minus.ptr = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if (ref_ft.ft == NULL) {
            if ((ref_ft_minus.name =
                 (char **) calloc(mLen,
                                  sizeof(*(ref_ft_minus.name)))) == NULL) {
                perror("process_sga: malloc");
                exit(1);
            }
        }
        /* Local Maxima Arrays  (position, count, feature name) */
        if ((lm_plus.cnt = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((lm_plus.pos =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((lm_minus.cnt = (int *) calloc(mLen, sizeof(int))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if ((lm_minus.pos =
             (unsigned long *) calloc(mLen,
                                      sizeof(unsigned long))) == NULL) {
            perror("process_sga: malloc");
            exit(1);
        }
        if (ref_ft.ft == NULL) {
            if ((lm_plus.name =
                 (char **) calloc(mLen,
                                  sizeof(*(lm_plus.name)))) == NULL) {
                perror("process_sga: malloc");
                exit(1);
            }
            if ((lm_minus.name =
                 (char **) calloc(mLen,
                                  sizeof(*(lm_minus.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 i = 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 (i >= SEQ_ID) {
                fprintf(stderr, "Seq ID is too long \"%s\" \n", buf);
                exit(1);
            }
            seq_id[i++] = *buf++;
        }
        while (isspace(*buf))
            buf++;
        /* FEATURE */
        i = 0;
        while (*buf != 0 && !isspace(*buf)) {
            if (i >= FT_MAX) {
                fprintf(stderr, "Feature is too long \"%s\" \n", buf);
                exit(1);
            }
            ft[i++] = *buf++;
        }
        while (isspace(*buf))
            buf++;
        /* Position */
        i = 0;
        while (isdigit(*buf)) {
            if (i >= POS_MAX) {
                fprintf(stderr, "Position is too large \"%s\" \n", buf);
                exit(1);
            }
            position[i++] = *buf++;
        }
        position[i] = 0;
        pos = atoi(position);
        while (isspace(*buf))
            buf++;
        /* Strand */
        strand = *buf++;
        while (isspace(*buf))
            buf++;
        /* Counts */
        i = 0;
        while (isdigit(*buf)) {
            if (i >= CNT_MAX) {
                fprintf(stderr, "Count is too large \"%s\" \n", buf);
                exit(1);
            }
            count[i++] = *buf++;
        }
        count[i] = 0;
        cnt = atoi(count);
        while (isspace(*buf))
            buf++;
        /* SGA Extension */
        i = 0;
        while (*buf != 0) {
            if (i >= EXT_MAX) {
                fprintf(stderr, "Extension is too long \"%s\" \n", buf);
                exit(1);
            }
            ext[i++] = *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 (!options.oriented) {
            if (k >= mLen - 1) {
                mLen *= 2;
#ifdef DEBUG
                fprintf(stderr,
                        "reallocating memory for ref_ft.pos ref_ft.name ref_ft.cnt (k=%d, size=%d)\n",
                        k, mLen);
#endif
                if ((ref_ft.pos =
                     (unsigned long *) realloc(ref_ft.pos,
                                               mLen *
                                               sizeof(unsigned long))) ==
                    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.ft == NULL) {
                    if ((ref_ft.name =
                         (char **) realloc(ref_ft.name,
                                           mLen *
                                           sizeof(*(ref_ft.name)))) ==
                        NULL) {
                        perror("process_sga: malloc");
                        exit(1);
                    }
                }
            }
        } else {
            if (k >= mLen - 1) {
                mLen *= 2;
#ifdef DEBUG
                fprintf(stderr,
                        "reallocating memory for ref_ft_plus.pos ref_ft_plus.name ref_ft_plus.cnt (k=%d, size=%d)\n",
                        k, mLen);
#endif
                if ((ref_ft_plus.pos =
                     (unsigned long *) realloc(ref_ft_plus.pos,
                                               mLen *
                                               sizeof(unsigned long))) ==
                    NULL) {
                    perror("process_sga: realloc");
                    exit(1);
                }
                if ((ref_ft_plus.cnt =
                     (int *) realloc(ref_ft_plus.cnt,
                                     mLen * sizeof(int))) == NULL) {
                    perror("process_sga: realloc");
                    exit(1);
                }
                if (ref_ft.ft == NULL) {
                    if ((ref_ft_plus.name =
                         (char **) realloc(ref_ft_plus.name,
                                           mLen *
                                           sizeof(*(ref_ft_plus.name)))) ==
                        NULL) {
                        perror("process_sga: malloc");
                        exit(1);
                    }
                }
            }
            if (j >= mLen1 - 1) {
                mLen1 *= 2;
#ifdef DEBUG
                fprintf(stderr,
                        "reallocating memory for ref_ft_minus.pos ref_ft_minus.name ref_ft_minus.cnt (j=%d, size=%d)\n",
                        j, mLen1);
#endif
                if ((ref_ft_minus.pos =
                     (unsigned long *) realloc(ref_ft_minus.pos,
                                               mLen1 *
                                               sizeof(unsigned long))) ==
                    NULL) {
                    perror("process_sga: realloc");
                    exit(1);
                }
                if ((ref_ft_minus.cnt =
                     (int *) realloc(ref_ft_minus.cnt,
                                     mLen1 * sizeof(int))) == NULL) {
                    perror("process_sga: realloc");
                    exit(1);
                }
                if (ref_ft.ft == NULL) {
                    if ((ref_ft_minus.name =
                         (char **) realloc(ref_ft_minus.name,
                                           mLen1 *
                                           sizeof(*(ref_ft_minus.name))))
                        == NULL) {
                        perror("process_sga: malloc");
                        exit(1);
                    }
                }
            }
        }
        /* Check Chromosome BEGINNING, process previous signal peaks and printout results */
        if (strcmp(seq_id, seq_id_prev) != 0) {
            if (!options.oriented) {
                ref_ft.pos[0] = ref_ft.pos[1] - Win1 / 2 - 1;
                ref_ft.pos[k + 1] = ref_ft.pos[k] + Win1 / 2 + 1;
                locate_peaks((int) k);
                if (ref_ft.ft == NULL) {
                    for (int i = 1; i <= (int) k; i++) {
                        free(ref_ft.name[i]);
                    }
                }
                k = 0;
            } else {
                /*printf("SEQID %s : k = %d j = %d calling locate_peaks_oriented\n", seq_id_prev, k, j); */
                ref_ft_plus.pos[0] = ref_ft_plus.pos[1] - Win1 / 2 - 1;
                ref_ft_plus.pos[k + 1] = ref_ft_plus.pos[k] + Win1 / 2 + 1;
                k1 = locate_peaks_oriented(&ref_ft_plus, &lm_plus, '+',
                                           (int) k);
                /*printf("after locate_peaks_oriented + : k1 = %d\n", k1); */
                ref_ft_minus.pos[0] = ref_ft_minus.pos[1] - Win1 / 2 - 1;
                ref_ft_minus.pos[j + 1] =
                    ref_ft_minus.pos[j] + Win1 / 2 + 1;
                j1 = locate_peaks_oriented(&ref_ft_minus, &lm_minus, '-',
                                           (int) j);
                /*printf("after locate_peaks_oriented - : j1 = %d\n", j1); */
                /*printf("calling merge_peaks: k1 %d, j1 %d, seq_id %s feat %s\n", k1, j1, seq_id_prev, ref_ft.ft); */
                merge_peaks((int) k1, (int) j1, seq_id_prev, ref_ft.ft);
                if (ref_ft.ft == NULL) {
                    for (int i = 1; i <= (int) k; i++) {
                        free(ref_ft_plus.name[i]);
                    }
                    for (int i = 1; i <= (int) j; i++) {
                        free(ref_ft_minus.name[i]);
                    }
                    for (int i = 0; i < (int) k1; i++) {
                        free(lm_plus.name[i]);
                    }
                    for (int i = 0; i < (int) j1; i++) {
                        free(lm_minus.name[i]);
                    }
                }
                k = 0;
                j = 0;
            }
            strcpy(seq_id_prev, seq_id);
        }
        if (!options.oriented) {
            if (ref_ft.ft == NULL) {
                k++;
                strcpy(ref_ft.seq_id, seq_id);
                ref_ft.name[k] = malloc(strlen(ft) + 3);
                strcpy(ref_ft.name[k], ft);
                strcat(ref_ft.name[k], "_p");
                ref_ft.pos[k] = (unsigned long) pos;
                if (cnt > Coff)
                    ref_ft.cnt[k] = Coff;
                else
                    ref_ft.cnt[k] = cnt;
            } else if (ref_ft.ft_str == '\0') {
                if (strcmp(ft, ref_ft.ft) == 0) {
                    k++;
                    strcpy(ref_ft.seq_id, seq_id);
                    /*strcpy(ref_ft.ft, ft); */
                    ref_ft.pos[k] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft.cnt[k] = Coff;
                    else
                        ref_ft.cnt[k] = cnt;
                }
            } else if (strand_flag == 1) {
                if (strand == ref_ft.ft_str) {
                    k++;
                    strcpy(ref_ft.seq_id, seq_id);
                    strcpy(ref_ft.ft, ft);
                    ref_ft.pos[k] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft.cnt[k] = Coff;
                    else
                        ref_ft.cnt[k] = cnt;
                }
            } else {
                if (strcmp(ft, ref_ft.ft) == 0 && strand == ref_ft.ft_str) {
                    k++;
                    strcpy(ref_ft.seq_id, seq_id);
                    /*strcpy(ref_ft.ft, ft); */
                    ref_ft.pos[k] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft.cnt[k] = Coff;
                    else
                        ref_ft.cnt[k] = cnt;
                }
            }
        } else {                /* Oriented */
            if (ref_ft.ft == NULL) {
                if (strand == '+') {
                    k++;
                    strcpy(ref_ft_plus.seq_id, seq_id);
                    ref_ft_plus.name[k] = malloc(strlen(ft) + 3);
                    strcpy(ref_ft_plus.name[k], ft);
                    strcat(ref_ft_plus.name[k], "_p");
                    ref_ft_plus.pos[k] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft_plus.cnt[k] = Coff;
                    else
                        ref_ft_plus.cnt[k] = cnt;
                } else if (strand == '-') {
                    j++;
                    strcpy(ref_ft_minus.seq_id, seq_id);
                    ref_ft_minus.name[j] = malloc(strlen(ft) + 3);
                    strcpy(ref_ft_minus.name[j], ft);
                    strcat(ref_ft_minus.name[j], "_p");
                    ref_ft_minus.pos[j] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft_minus.cnt[j] = Coff;
                    else
                        ref_ft_minus.cnt[j] = cnt;
                }
            } else {
                if (strcmp(ft, ref_ft.ft) == 0 && strand == '+') {
                    k++;
                    strcpy(ref_ft_plus.seq_id, seq_id);
                    ref_ft_plus.pos[k] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft_plus.cnt[k] = Coff;
                    else
                        ref_ft_plus.cnt[k] = cnt;
                } else if (strcmp(ft, ref_ft.ft) == 0 && strand == '-') {
                    j++;
                    strcpy(ref_ft_minus.seq_id, seq_id);
                    ref_ft_minus.pos[j] = (unsigned long) pos;
                    if (cnt > Coff)
                        ref_ft_minus.cnt[j] = Coff;
                    else
                        ref_ft_minus.cnt[j] = cnt;
                }
            }
        }
    }                           /* End of While */
    free(s);
    /* Locate signal peaks for the last chromosome */
    if (!options.oriented) {
        ref_ft.pos[0] = ref_ft.pos[1] - Win1 / 2 - 1;
        ref_ft.pos[k + 1] = ref_ft.pos[k] + Win1 / 2 + 1;
        locate_peaks((int) k);
        if (ref_ft.ft == NULL) {
            for (int i = 1; i <= (int) k; i++) {
                free(ref_ft.name[i]);
            }
        }
        free(ref_ft.name);
        free(ref_ft.pos);
        free(ref_ft.cnt);
        free(ref_ft.nct);
        free(ref_ft.npo);
        free(ref_ft.ptr);

        free(lm_plus.pos);
        free(lm_plus.cnt);
        free(lm_minus.pos);
        free(lm_minus.cnt);
    } else {
        ref_ft_plus.pos[0] = ref_ft_plus.pos[1] - Win1 / 2 - 1;
        ref_ft_plus.pos[k + 1] = ref_ft_plus.pos[k] + Win1 / 2 + 1;
        k1 = locate_peaks_oriented(&ref_ft_plus, &lm_plus, '+', (int) k);
        ref_ft_minus.pos[0] = ref_ft_minus.pos[1] - Win1 / 2 - 1;
        ref_ft_minus.pos[j + 1] = ref_ft_minus.pos[j] + Win1 / 2 + 1;
        j1 = locate_peaks_oriented(&ref_ft_minus, &lm_minus, '-', (int) j);
        merge_peaks((int) k1, (int) j1, seq_id_prev, ref_ft.ft);
        if (ref_ft.ft == NULL) {
            for (int i = 1; i <= (int) k; i++) {
                free(ref_ft_plus.name[i]);
            }
            for (int i = 1; i <= (int) j; i++) {
                free(ref_ft_minus.name[i]);
            }
            for (int i = 0; i < (int) k1; i++) {
                free(lm_plus.name[i]);
            }
            for (int i = 0; i < (int) j1; i++) {
                free(lm_minus.name[i]);
            }
        }
        free(ref_ft_plus.name);
        free(ref_ft_minus.name);

        free(ref_ft_plus.pos);
        free(ref_ft_plus.cnt);
        free(ref_ft_minus.pos);
        free(ref_ft_minus.cnt);

        free(ref_ft_plus.nct);
        free(ref_ft_plus.npo);
        free(ref_ft_plus.ptr);
        free(ref_ft_minus.nct);
        free(ref_ft_minus.npo);
        free(ref_ft_minus.ptr);

        free(lm_plus.name);
        free(lm_plus.pos);
        free(lm_plus.cnt);
        free(lm_minus.name);
        free(lm_minus.pos);
        free(lm_minus.cnt);
    }
    if (input != stdin) {
        fclose(input);
    }
    return 0;
}

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

    while (1) {
        int c = getopt(argc, argv, "f:di:horw:v:t:c:");
        if (c == -1)
            break;
        switch (c) {
        case 'c':
            Coff = atoi(optarg);
            break;
        case 'd':
            options.debug = 1;
            break;
        case 'f':
            Feature = optarg;
            break;
        case 'h':
            options.help = 1;
            break;
        case 'i':
            options.dbPath = optarg;
            options.db = 1;
            break;
        case 'o':
            options.oriented = 1;
            break;
        case 'r':
            options.refine = 1;
            break;
        case 'w':
            Win1 = atoi(optarg);
            break;
        case 'v':
            Win2 = atoi(optarg);
            break;
        case 't':
            Thres = atoi(optarg);
            break;
        default:
            printf("?? getopt returned character code 0%o ??\n", c);
        }
    }
    if (optind > argc || options.help == 1 || Win1 == 0 || Win2 == 0
        || Coff < 0) {
        show_help();
        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);
                show_help();
                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, " Integration range (Window) : %d\n\n", Win1);
        fprintf(stderr, " Minimal distance (Vicinity) : %d\n\n", Win2);
        fprintf(stderr, " Peak Threshold : %d\n\n", Thres);
    }
    /* 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 + 2) * 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';
        if (!options.oriented) {
            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 (!options.oriented) {
        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_size() == 0) {
        if (options.debug)
            fprintf(stderr,
                    " HASH Table for chromosome size initialized\n");
    } else {
        return 1;
    }
    if (process_sga(input, argv[optind++]) != 0) {
        return 1;
    }
    if (ref_ft.ft != NULL) {
        free(ref_ft.ft);
    }
    return 0;
}
