/*
  sga2wig.c

  Convert SGA file to Wiggle Track (WIG) format.

  There are two options for formatting wiggle data:
    - variableStep [0]
    - fixedStep    [1]

   variableStep is for data with irregular intervals between new data points.
   It begins with a declaration line and is followed by two columns containing chromosome positions and data values:

       variableStep  chrom=chrN  [span=windowSize]
       chromStartA  dataValueA
       chromStartB  dataValueB
       ... etc ...  ... etc ...

  fixedStep is for data with regular intervals between new data values.
  It begins with a declaration line and is followed by a single column of data values:

         fixedStep  chrom=chrN  start=position  step=stepInterval  [span=windowSize = stepInterval]
         dataValue1
         dataValue2
         ... etc ...
  (For more details on WIG format see https://genome.ucsc.edu/goldenPath/help/wiggle.html)

  # Arguments:
  #   SGA File
  # Options:
  #   Set Path to locate chr_NC_gi file
  #   (used for NCBI id to chromosome name conversion)
  #   Set WIG format option: variableStep[0]/fixedStep[1] (def=0)
  #   Set Chromosome number  (by def=0 i.e. take all chromosomes))
  #   Set chromosome start coordinate (def=-1, i.e. entire chromosome region)
  #   Set chromosome end coordinate (def=-1, i.e. entire chromosome region)
  #   Set count cut-off for SGA count field (def=99999)
  #   Set normalisation factor for total raead counts within a step interval
  #    - For fixedStep WIG format (def=0, i.e. take sum of total counts)
  #   Set Wiggle span parameter (def=150)
  #   Set Wiggle step parameter (fixedStep Format) (def=150 step=span)
  #   SGA input file is a peak file, expand peak coordinates by span value
  #   on both upstream and downstream directions
  #   Set data viewing paramenter: autoscale (def=OFF)
  #   Set data viewing paramenter: always0 (include zero) (def=OFF)
  #   Set data viewing paramenter: wfunction function (def=mean+whiskers, maximum|mean|minimum)
  #   Set data viewing paramenter: smoothing window (def=OFF[0], <grade>=0,2..16)
  #   Set data viewing paramenter: visibility mode (def=full, dense|hide)
  #   Set Wiggle track name (def=Custom-Wig)
  #   Set Wiggle track description (def=ChIP-Seq)
  #   Set Wiggle track color (def= 0,200,100)

  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 <getopt.h>
#include <assert.h>
#include <math.h>
#include <sys/stat.h>
#include "hashtable.h"
#ifdef DEBUG
#include <mcheck.h>
#endif

#include "version.h"

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

typedef struct _options_t {
    char *dbPath;
    int help;
    int debug;
    int db;
    int wigFormat;
    char *chrNb;
    char *chrName;
    int chrFlag;
    unsigned long chrStart;
    unsigned long chrEnd;
    int cutOff;
    int norm;
    int span;
    int step;
    int peakFlag;
    char *autoscale;
    char *always0;
    char *wfunction;
    char *smoothing;
    char *visibility;
    char *trackName;
    char *trackDesc;
    char *trackColor;
} options_t;

static options_t options;

typedef struct _feature_t {
    char seq_id[SEQ_ID];
    unsigned long *pos;
    int *cnt;
} feature_t, *feature_p_t;

static feature_t sga_ft;

static hash_table_t *ac_table = NULL;
static hash_table_t *size_table = NULL;

unsigned long endPos = 0;
unsigned long startPos = 0;

char *progname = NULL;

void show_help(void) {
    fprintf(stderr, "Usage: %s [options] [<] <SGA file|stdin>\n"
            "      - version %s\n"
            "      where options are:\n"
            "  \t\t -d|--debug                Produce Debug information\n"
            "  \t\t -h|--help                 Show this Help text\n"
            "  \t\t -i|--db <path>            Use <path> to locate the chr_NC_gi and chr_size files\n"
            "  \t\t                           [default is: $CWD]\n"
            "  \t\t --format <0|1>            Set Wiggle Track data format: variableStep[def=0]/fixedStep[1]\n"
            "  \t\t --chr       <int>         Chromosome number [def: 0 (all chromosomes)]\n"
            "  \t\t --chrnb     <int>         Deprecated, please use --chr\n"
            "  \t\t -b|--start  <int>         Chromosome start [def: -1 (entire chromosome)]\n"
            "  \t\t -e|--end    <int>         Chromosome end [def: -1 (entire chromosome)]\n"
            "  \t\t -C|--coff   <int>         Count cut-off for the SGA input file [def=99999]\n"
            "  \t\t -s|--span   <int>         Wiggle Track Span(/stepInterval) parameter [def=%d]\n"
            "  \t\t                           For fixedStep data format, it defines the step parameter\n"
            "  \t\t -n|--norm   <int>         Normalization factor for total tag counts within step intervals [def=0]\n"
            "  \t\t                           This option is only valid for fixedStep data format\n"
            "  \t\t --peakf                   Indicate that the The SGA input file represents a peak file\n"
            "  \t\t                           [i.e. coordinates are peak centers]\n"
            "  \t\t                           In such case, the span range begins upstream of [span=]%d bp\n"
            "  \t\t                           chromosome position specified, and ends [span=]%d bp downstream\n"
            "  \t\t --name <name>             Set name for track name field [def. name=SGA-feature]\n"
            "  \t\t --desc <desc>             Set track description field [def. desc=\"ChIP-Seq Custom data\"]\n"
            "  \t\t --color <col>             Define the track color in comma-separated RGB values [def. 100,100,100]\n"
            "  \t\t --autoscale  <on|off>     Data viewing paramenter: set auto-scale to UCSC data view [def=OFF]\n"
            "  \t\t --always0    <on|off>     Data viewing paramenter: always include zero [def=OFF]\n"
            "  \t\t --wfunction  <func>       Data viewing paramenter: windowing function [def=mean+whiskers|maximum|mean|minimum]\n"
            "  \t\t --smoothing  <grade>      Data viewing paramenter: smoothing window [def=OFF[0], <grade>=0,2..16]\n"
            "  \t\t --visibility <mode>       Display mode: [def=full|dense|hide]\n"
            "\n\tConvert SGA format into Wiggle Track format (WIG).\n"
            "\tWIG format is line-oriented, and is composed of declaration lines and data lines.\n\n"
            "\t- variableStep is for data with irregular intervals between new data points.\n"
            "\t  It begins with a declaration line and is followed by two columns containing chromosome positions and data values:\n\n"
            "\t  variableStep  chrom=chrN  [span=windowSize]\n"
            "\t  chromStartA  dataValueA\n"
            "\t  chromStartB  dataValueB\n"
            "\t  ... etc ...  ... etc ...\n\n"
            "\t- fixedStep is for data with regular intervals between new data values.\n"
            "\t  It begins with a declaration line and is followed by a single column of data values:\n\n"
            "\t  fixedStep  chrom=chrN  start=position  step=stepInterval  [span=windowSize = stepInterval]\n"
            "\t  dataValue1\n"
            "\t  dataValue2\n"
            "\t  ... etc ...\n\n",
            progname, VERSION, options.span, options.span, options.span);
}

int process_ac(void) {
    FILE *input;
    int c;
    char buf[LINE_SIZE];
    char *chrFile;
    char chrom[18];
    int cLen;

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

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

    do {
        c = fgetc(input);
    } while (c != '\n');

    ac_table = hash_table_new(MODE_COPY);
    while (fgets(buf, LINE_SIZE, input) != NULL) {
        char *s;
        char chr_nb[CHR_NB] = "";
        char ncbi_ac[AC_MAX] = "";
        int i = 0;
        int nb_len = 0;
        int ac_len = 0;
        /*int valid = 1; */
        s = buf;
        /* Check line */
        /* Get first character: if # skip line */
        if (*s == '#')
            continue;
        /* Chrom NB */
        while (*s != 0 && !isspace(*s)) {
            if (i >= CHR_NB) {
                fprintf(stderr, "AC too long in %s\n", s);
                fclose(input);
                exit(1);
            }
            chr_nb[i++] = *s++;
        }
        if (i < CHR_NB)
            chr_nb[i] = 0;
        nb_len = i + 1;
        while (isspace(*s))
            s++;
        /* 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;
        strcpy(chrom, "chr");
        strcat(chrom, chr_nb);
        nb_len = (int) strlen(chrom) + 1;
        /* printf("adding key %s (len = %d) value %s (ac) (len = %d) to hash table\n", chrom, nb_len, ncbi_ac, ac_len); */
        /* Store both NCBI identifier to chrom number and chrom number to chrom number keys */
        hash_table_add(ac_table, ncbi_ac, (size_t) ac_len, chrom,
                       (size_t) nb_len);
        if (options.debug) {
            char *cn =
                hash_table_lookup(ac_table, ncbi_ac, (size_t) ac_len);
            fprintf(stderr,
                    " AC Hash table: %s (len = %d) -> %s (len = %d)\n",
                    ncbi_ac, ac_len, cn, nb_len);
        }
        hash_table_add(ac_table, chrom, (size_t) nb_len, chrom,
                       (size_t) nb_len);
        if (options.debug) {
            char *cn = hash_table_lookup(ac_table, chrom, (size_t) nb_len);
            fprintf(stderr,
                    " AC Hash table: %s (len = %d) -> %s (len = %d)\n",
                    chrom, nb_len, cn, nb_len);
        }
    }
    return 0;
}

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 = 21 + 10;
        if ((chrSizeFile = (char *) malloc(cLen * sizeof(char))) == NULL) {
            perror("process_ac: malloc");
            exit(1);
        }
        /*
         * Replace /home/local/db/genome (which cannot be assumed to exist on users'
         * machines) with $CWD.
         strcpy(chrSizeFile, "/home/local/db/genome");
         */
        strcpy(chrSizeFile, ".");
    }
    strcat(chrSizeFile, "/chr_size");

    input = fopen(chrSizeFile, "r");
    if (input == NULL) {
        fprintf(stderr, "Could not open file %s: %s(%d)\n",
                chrSizeFile, strerror(errno), errno);
        return 1;
    }

    do {
        c = fgetc(input);
    } while (c != '\n');

    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
write_wig_vs(unsigned long *pos, int *counts, unsigned int len, char *ac) {
    unsigned int i = 0;
    unsigned int k = 0;
    unsigned long *p = pos;
    int *cnt = counts;
    int ac_len = 0;

    if (endPos == 0) {
        return;
    }
    if (ac == NULL)
        return;
    ac_len = (int) strlen(ac) + 1;
    char *chr_nb = hash_table_lookup(ac_table, ac, (size_t) ac_len);
    if (chr_nb == NULL)
        return;
    printf("variableStep chrom=%s span=%d\n", chr_nb, options.span);
    printf("1 0\n");
    int offset = 2;
    for (i = 0; i < len; i++) {
        if (options.peakFlag) {
            long new_pos = *p++ - options.span / 2;
            if (new_pos <= 1) {
                new_pos = offset++;
            }
            printf("%ld %d\n", new_pos, *cnt++);
        } else {
            k = 1;
            /* allows data composed of contiguous runs of options.span bases
               the same data value to be specified more succinctly */
            while ((i + k < len) && (p[i + k] < p[i] + options.span)) {
                if (cnt[i] == cnt[i + k])
                    k++;
                else
                    break;
            }
            //printf("%lu %d\n", *s++, *cnt++);
            printf("%lu %d\n", p[i], cnt[i]);
            if (k == (unsigned) options.span)
                i = i + k - 1;
        }
    }
}

void
write_wig_fs(unsigned long *pos, int *counts, unsigned int len, char *ac) {
    unsigned int i = 0;
    unsigned int k = 0;
    //unsigned long *s = pos;
    //int *cnt = counts;
    int ac_len = 0;
    char *csize = NULL;
    unsigned long chr_size = 0;
    long sum = 0;

    if (endPos == 0) {
        return;
    }
    unsigned long end_loop = endPos - options.step;
    if (ac == NULL)
        return;
    ac_len = (int) strlen(ac) + 1;
    csize = hash_table_lookup(size_table, ac, (size_t) ac_len);
    if (csize == NULL)
        return;
    chr_size = (unsigned long) atoi(csize);
    char *chr_nb = hash_table_lookup(ac_table, ac, (size_t) ac_len);
    if (chr_nb == NULL)
        return;
    if (end_loop > chr_size) {
        end_loop = chr_size - options.span + 1;
    }
    printf("fixedStep chrom=%s start=%lu step=%d span=%d\n", chr_nb,
           startPos, options.step, options.span);
    for (i = startPos; i <= end_loop; i += options.step) {
        while ((k < len) && (pos[k] <= i + options.step)) {
            sum += counts[k++];
        }
        if (options.norm) {
            sum = (long) ((sum * options.norm) / options.step);
        }
        printf("%ld\n", sum);
        sum = 0;
    }
}

int process_sga(FILE *input, char *iFile) {
    //unsigned long start, end;
    unsigned long pos;
    int cnt;
    //char *name;
    int first = 1;
    int wig_hdr = 1;
    char *s, *res, *buf;
    size_t bLen = LINE_SIZE;
    size_t sga_mLen = BUF_SIZE;
    char seq_id[SEQ_ID] = "";
    char seq_id_prev[SEQ_ID] = "";
    unsigned int k = 0;

    if (options.debug)
        fprintf(stderr, " Processing SGA file %s\n", iFile);
    if ((s = malloc(bLen * sizeof(char))) == NULL) {
        perror("process_sga: malloc");
        exit(1);
    }

    /* Malloc position and count arrays */
    if ((sga_ft.pos =
         (unsigned long *) calloc(sga_mLen,
                                  sizeof(unsigned long))) == NULL) {
        perror("process_sga: malloc");
        exit(1);
    }
    if ((sga_ft.cnt = (int *) calloc(sga_mLen, sizeof(int))) == NULL) {
        perror("process_sga: malloc");
        exit(1);
    }
#ifdef DEBUG
    int c = 1;
#endif
    while ((res = fgets(s, (int) bLen, input)) != NULL) {
        size_t cLen = strlen(s);
        char ft[FT_MAX] = "";
        char position[POS_MAX] = "";
        char count[CNT_MAX] = "";
#ifdef DEBUG
        char strand = '\0';
#endif
        char ext[EXT_MAX];
        int i = 0;
        int id_len = 0;
        char *cn = NULL;

        memset(ext, 0, (size_t) EXT_MAX);
        memset(seq_id, 0, (size_t) SEQ_ID);
        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 = *buf++;
        buf++;
#ifdef DEBUG
        /* Strand */
        strand = *buf;
#endif
        while (isspace(*buf))
            buf++;
        /* Counts */
        i = 0;
        while (isdigit(*buf) || *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
        //printf(" seq ID: %s   Feat: %s (%c)  Pos: %d  Cnts: %d Ext: %s\n", seq_id, ft, strand, pos, cnt, ext);
        if (k > sga_mLen - 1) {
            sga_mLen *= 2;
            if ((sga_ft.pos =
                 (unsigned long *) realloc(sga_ft.pos,
                                           sga_mLen *
                                           sizeof(unsigned long))) ==
                NULL) {
                perror("process_sga: realloc pos array");
                exit(1);
            }
            if ((sga_ft.cnt =
                 (int *) realloc(sga_ft.cnt,
                                 sga_mLen * sizeof(int))) == NULL) {
                perror("process_sga: realloc count array");
                exit(1);
            }
        }
        if (strcmp(seq_id, seq_id_prev) != 0) {
            if (k > 0)
                endPos = sga_ft.pos[k - 1];
            else
                endPos = 0;
            if (options.wigFormat) {
                write_wig_fs(sga_ft.pos, sga_ft.cnt, k, seq_id_prev);
            } else {
                write_wig_vs(sga_ft.pos, sga_ft.cnt, k, seq_id_prev);
            }
            strcpy(seq_id_prev, seq_id);
            k = 0;
            wig_hdr = 1;
        }
        /* Set WIG annotation track line */
        if (first) {
            if (options.trackName == NULL) {
                if ((options.trackName =
                     malloc((strlen(ft) + 1) * sizeof(char))) == NULL) {
                    perror("process_sga: malloc Feature");
                    exit(1);
                }
                strcpy(options.trackName, ft);
            }
            if (options.chrName != NULL) {
                //printf("chr : %s start : %d  end : %d\n", options.chrName, options.chrStart, options.chrEnd);
                if ((options.chrStart != -1lu) && (options.chrEnd != -1lu)) {
                    printf("browser position %s:%lu-%lu\n",
                           options.chrName, options.chrStart,
                           options.chrEnd);
                } else if (options.chrStart != -1lu) {
                    printf("browser position %s:%lu-%lu\n",
                           options.chrName, options.chrStart,
                           options.chrStart + 100000);
                } else {
                    printf("browser position %s:%d-%d\n", options.chrName,
                           1, 100000);
                }
            } else {
                printf("browser full refGene\n");
            }
            printf
                ("track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=%s color=%s autoScale=%s alwaysZero=%s maxHeightPixels=100:50:20 graphType=bar priority=30 windowingFunction=%s smoothingWindow=%s\n",
                 options.trackName, options.trackDesc, options.visibility,
                 options.trackColor, options.autoscale, options.always0,
                 options.wfunction, options.smoothing);
            first = 0;
        }
        /* Get Chromosome name */
        //printf("Seq ID: %s pos: %lu  cnt: %d\n", seq_id, pos, cnt);
        id_len = (int) strlen(seq_id) + 1;
        cn = hash_table_lookup(ac_table, seq_id, (size_t) id_len);
        //printf("Chr name: %s\n", cn);
        if (options.chrFlag) {  /* Chromosome name has been specified */
            if ((options.chrStart != -1lu) && (options.chrEnd != -1lu)) {
                if (cn != NULL) {
                    //printf ("chr: %s\tselected chr: %s\n", cn, options.chrName);
                    if (strcmp(cn, options.chrName) == 0
                        && pos >= options.chrStart
                        && pos <= options.chrEnd) {
                        //printf("Adding chr: %s  pos: %lu ([%lu-%lu])\n", cn, pos, options.chrStart, options.chrEnd);
                        if (cnt > options.cutOff)
                            sga_ft.cnt[k] = options.cutOff;
                        else
                            sga_ft.cnt[k] = cnt;
                        sga_ft.pos[k] = pos;
                        k++;
                        if (wig_hdr) {
                            startPos = options.chrStart;
                            wig_hdr = 0;
                        }
                    }
                }
            } else if (options.chrStart != -1lu) {
                if (cn != NULL) {
                    if (strcmp(cn, options.chrName) == 0
                        && pos >= options.chrStart) {
                        if (cnt > options.cutOff)
                            sga_ft.cnt[k] = options.cutOff;
                        else
                            sga_ft.cnt[k] = cnt;
                        sga_ft.pos[k] = pos;
                        k++;
                        if (wig_hdr) {
                            startPos = options.chrStart;
                            wig_hdr = 0;
                        }
                    }
                }
            } else {
                if (cn != NULL) {
                    if (strcmp(cn, options.chrName) == 0) {
                        if (cnt > options.cutOff)
                            sga_ft.cnt[k] = options.cutOff;
                        else
                            sga_ft.cnt[k] = cnt;
                        sga_ft.pos[k] = pos;
                        k++;
                        if (wig_hdr) {
                            startPos = pos;
                            wig_hdr = 0;
                        }
                    }
                }
            }
        } else {                /* Scan entire SGA */
            if (cnt > options.cutOff)
                sga_ft.cnt[k] = options.cutOff;
            else
                sga_ft.cnt[k] = cnt;
            sga_ft.pos[k] = pos;
            k++;
            if (wig_hdr) {
                startPos = pos;
                wig_hdr = 0;
            }
        }
    }                           /* End of While */
    if (input != stdin) {
        fclose(input);
    }
    /* The last time (at EOF) */
    if (k > 0)
        endPos = sga_ft.pos[k - 1];
    else
        return 0;
    //printf("The last time: SEQ_ID %s, endPos: %lu len: %d\n", seq_id, endPos, k);
    if (options.wigFormat) {
        write_wig_fs(sga_ft.pos, sga_ft.cnt, k, seq_id);
    } else {
        write_wig_vs(sga_ft.pos, sga_ft.cnt, k, seq_id);
    }
    free(sga_ft.pos);
    free(sga_ft.cnt);
    return 0;
}

int main(int argc, char *argv[]) {
    progname = argv[0];
    options.wigFormat = 0;
    options.peakFlag = 0;
    options.chrNb = NULL;
    options.chrName = NULL;
    options.chrFlag = 0;
    options.chrStart = -1;
    options.chrEnd = -1;
    options.cutOff = 99999;
    options.span = 150;
    options.step = 150;
#ifdef DEBUG
    mcheck(NULL);
    mtrace();
#endif
    FILE *input;
    //int i = 0;
    //int j = 0;
    int usage = 0;

    static struct option long_options[] = {
        /* These options may or may not set a flag.
           We distinguish them by their indices. */
        { "debug", no_argument, 0, 'd' },
        { "help", no_argument, 0, 'h' },
        { "db", required_argument, 0, 'i' },
        { "format", required_argument, 0, 0 },
        { "chrnb", required_argument, 0, 0 },
        { "chr", required_argument, 0, 0 },
        { "start", required_argument, 0, 'b' },
        { "end", required_argument, 0, 'e' },
        { "coff", required_argument, 0, 'C' },
        { "norm", required_argument, 0, 'n' },
        { "span", required_argument, 0, 's' },
        { "name", required_argument, 0, 0 },
        { "desc", required_argument, 0, 0 },
        { "color", required_argument, 0, 0 },
        { "autoscale", required_argument, 0, 0 },
        { "always0", required_argument, 0, 0 },
        { "wfunction", required_argument, 0, 0 },
        { "smoothing", required_argument, 0, 0 },
        { "visibility", required_argument, 0, 0 },
        /* These option only sets a flag. */
        { "peakf", no_argument, &options.peakFlag, 1 },
        { 0, 0, 0, 0 }
    };

    int option_index = 0;

    while (1) {
        int c = getopt_long(argc, argv, "dhi:o:n:b:e:C:f:s:", long_options,
                            &option_index);
        if (c == -1)
            break;
        switch (c) {
        case 'd':
            options.debug = 1;
            break;
        case 'h':
            options.help = 1;
            break;
        case 'i':
            options.dbPath = optarg;
            options.db = 1;
            break;
        case 'o':
            fprintf(stderr,
                    "WARNING: option -o is deprecated; please use --format instead.\n");
            options.wigFormat = atoi(optarg);
            break;
        case 'b':
            options.chrStart = atoi(optarg);
            break;
        case 'e':
            options.chrEnd = atoi(optarg);
            if (options.chrEnd != -1lu && options.chrStart == -1lu)
                options.chrStart = 1;
            break;
        case 'C':
            options.cutOff = atoi(optarg);
            break;
        case 'n':
            options.norm = atoi(optarg);
            break;
        case 's':
            options.span = atoi(optarg);
            options.step = options.span;
            break;
        case 0:
            /* This option is to set the annotation track line */
            if (strcmp(long_options[option_index].name, "name") == 0) {
                options.trackName = optarg;
            }
            if (strcmp(long_options[option_index].name, "chrnb") == 0
                || strcmp(long_options[option_index].name, "chr") == 0) {
                if (strcmp(long_options[option_index].name, "chrnb") == 0)
                    fprintf(stderr,
                            "WARNING: option --chrnb is deprecated; please use --chr.\n");
                options.chrNb = optarg;
                if (strcmp(options.chrNb, "0") != 0) {
                    if ((options.chrName =
                         malloc((strlen(options.chrNb) +
                                 4) * sizeof(char))) == NULL) {
                        perror("process_sga: malloc Chr Name");
                        exit(1);
                    }
                    strcpy(options.chrName, "chr");
                    strcat(options.chrName, options.chrNb);
                    options.chrFlag = 1;
                }
            }
            if (strcmp(long_options[option_index].name, "desc") == 0) {
                options.trackDesc = optarg;
            }
            if (strcmp(long_options[option_index].name, "color") == 0) {
                options.trackColor = optarg;
            }
            if (strcmp(long_options[option_index].name, "autoscale") == 0) {
                options.autoscale = optarg;
            }
            if (strcmp(long_options[option_index].name, "always0") == 0) {
                options.always0 = optarg;
            }
            if (strcmp(long_options[option_index].name, "wfunction") == 0) {
                options.wfunction = optarg;
            }
            if (strcmp(long_options[option_index].name, "smoothing") == 0) {
                options.smoothing = optarg;
            }
            if (strcmp(long_options[option_index].name, "visibility") == 0) {
                options.visibility = optarg;
            }
            if (strcmp(long_options[option_index].name, "format") == 0) {
                options.wigFormat = atoi(optarg);
            }
            break;
        default:
            printf("?? getopt returned character code 0%o ??\n", c);
            usage = 1;
        }
    }
    /* printf("optind: %d argc: %d\n", optind, argc); */
    if (optind > argc || options.help == 1 || usage) {
        show_help();
        return 1;
    }
    /* printf("argc: %d optind: %d\n", argc, optind); */
    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.trackDesc == NULL) {
        if ((options.trackDesc = malloc(32 * sizeof(char))) == NULL) {
            perror("main: malloc trackDesc");
            exit(1);
        }
        strcpy(options.trackDesc, "ChIP-Seq Custom data");
    }
    if (options.trackColor == NULL) {
        if ((options.trackColor = malloc(12 * sizeof(char))) == NULL) {
            perror("main: malloc trackColor");
            exit(1);
        }
        strcpy(options.trackColor, "0,200,100");
    }
    if (options.autoscale == NULL) {
        if ((options.autoscale = malloc(4 * sizeof(char))) == NULL) {
            perror("main: malloc autoscale");
            exit(1);
        }
        strcpy(options.autoscale, "off");
    }
    if (options.always0 == NULL) {
        if ((options.always0 = malloc(4 * sizeof(char))) == NULL) {
            perror("main: malloc always0");
            exit(1);
        }
        strcpy(options.always0, "off");
    }
    if (options.wfunction == NULL) {
        if ((options.wfunction = malloc(5 * sizeof(char))) == NULL) {
            perror("main: malloc wfunction");
            exit(1);
        }
        strcpy(options.wfunction, "mean");
    }
    if (options.smoothing == NULL) {
        if ((options.smoothing = malloc(4 * sizeof(char))) == NULL) {
            perror("main: malloc smoothing");
            exit(1);
        }
        strcpy(options.smoothing, "off");
    }
    if (options.visibility == NULL) {
        if ((options.visibility = malloc(5 * sizeof(char))) == NULL) {
            perror("main: malloc visibility");
            exit(1);
        }
        strcpy(options.visibility, "full");
    }
    if (options.peakFlag) {
        options.span *= 2;
    }
    if (options.debug) {
        fprintf(stderr, " Arguments:\n");
        fprintf(stderr, " Wiggle data format : %d\n", options.wigFormat);
        fprintf(stderr, " Chromosome Nb : %s  Name : %s\n", options.chrNb,
                options.chrName);
        fprintf(stderr, " Chromosome Start : %lu\n", options.chrStart);
        fprintf(stderr, " Chromosome End : %lu\n", options.chrEnd);
        fprintf(stderr, " Count cut-off : %d\n", options.cutOff);
        fprintf(stderr, " Normalisation factor for fixedStep WIG : %d\n",
                options.norm);
        fprintf(stderr, " WIG Span/Step Interval: %d\n", options.span);
        fprintf(stderr, " Peak Flag (for SGA input file): %d\n",
                options.peakFlag);
        fprintf(stderr, " Wiggle Track Display options:\n");
        fprintf(stderr, " -----------------------------\n");
        fprintf(stderr, " Track name: %s\n", options.trackName);
        fprintf(stderr, " Track description: %s\n", options.trackDesc);
        fprintf(stderr, " Track color: %s\n", options.trackColor);
        fprintf(stderr, " Autoscale: %s\n", options.autoscale);
        fprintf(stderr, " Always0: %s\n", options.always0);
        fprintf(stderr, " Wfunction: %s\n", options.wfunction);
        fprintf(stderr, " Smoothing: %s\n", options.smoothing);
        fprintf(stderr, " Visibility: %s\n", options.visibility);
        fprintf(stderr, "\n");
    }
    // See sga2bed.c, ll. 808-834
    // if (!((optind == 1) && (optind == argc))) {
    if (process_ac() == 0) {
        if (options.debug)
            fprintf(stderr,
                    " HASH Table for chromosome access identifier initialized\n");
    } else {
        return 1;
    }
    if (process_size() == 0) {
        if (options.debug)
            fprintf(stderr,
                    " HASH Table for chromosome size initialized\n");
    } else {
        return 1;
    }
    /* }
       else {
       return 1;
       }
     */
    if (process_sga(input, argv[optind++]) != 0) {
        return 1;
    }
    return 0;
}
