/*
  countsga.c

  Count tags and compute total sequence length of an SGA file

  # Arguments:
  # feature type, count cut-off
  # SGA file

  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>
#ifdef DEBUG
#include <mcheck.h>
#endif

#include "version.h"

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

typedef struct _options_t {
    int help;
    int debug;
    int cutOff;
    char *ftName;
    char strand;
} options_t;

static options_t options;
char *Feature = NULL;
int ft_specs = 1;

unsigned long Counts = 0;       /* Total Target Counts    */
unsigned long Len = 0;          /* Total Sequence Length  */

int process_sga(FILE *input) {
    char seq_id_prev[SEQ_ID] = "";
    unsigned long pos;
    int cnt = 0;
    unsigned long last_pos = 0;
    char *s, *res, *buf;
    size_t bLen = LINE_SIZE;

    if ((s = malloc(bLen * sizeof(char))) == NULL) {
        perror("process_sga: malloc");
        exit(1);
    }
    /*
       if (options.debug)
       fprintf(stderr, "Processing file %s\n", iFile);
     */
#ifdef DEBUG
    int lc = 1;
#endif
    /*
       while ((res = fgets(s, (int) bLen, f)) != NULL) {
     */
    while ((res = fgets(s, (int) bLen, input)) != NULL) {
        char seq_id[SEQ_ID] = "";
        char feature[FT_MAX] = "";
        char position[POS_MAX] = "";
        char count[CNT_MAX] = "";
        char strand = '\0';
        char ext[EXT_MAX];
        size_t cLen = strlen(s);
        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), f);
             */
            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);
            }
            feature[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 = (unsigned long) atol(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: %lu  Cnts: %d Ext: %s\n",
             lc++, seq_id, feature, strand, pos, cnt, ext);
#endif
        if (!strcmp(feature, "END")) {
            continue;
        }
        if (strcmp(seq_id, seq_id_prev) != 0) {
            Len += last_pos;
            strcpy(seq_id_prev, seq_id);
        }
        if (ft_specs == 0) {
            if (cnt > options.cutOff)
                Counts += (unsigned int) options.cutOff;
            else
                Counts += (unsigned int) cnt;
        } else {
            if (options.strand == '\0') {
                if (strcmp(feature, options.ftName) == 0) {
                    if (cnt > options.cutOff)
                        Counts += (unsigned int) options.cutOff;
                    else
                        Counts += (unsigned int) cnt;
                }
            } else {
                if (strcmp(feature, options.ftName) == 0
                    && strand == options.strand) {
                    if (cnt > options.cutOff)
                        Counts += (unsigned int) options.cutOff;
                    else
                        Counts += (unsigned int) cnt;
                }
            }
        }
        last_pos = pos;
    }                           /* End of While */
    free(s);
    /* The last time (at EOF) */
    Len += last_pos;
    /*
       fprintf (stderr, "Total Tag Counts : %lu , Total Sequence Len : %lu\n", Counts, Len);
     */
    printf("Total Tag Counts : %lu , Total Sequence Len : %lu\n", Counts,
           Len);
    return 0;
}

int main(int argc, char *argv[]) {
#ifdef DEBUG
    mcheck(NULL);
    mtrace();
#endif
    FILE *input;
    options.cutOff = 1;
    options.ftName = NULL;
    options.strand = '\0';

    while (1) {
        int c = getopt(argc, argv, "f:dhc:");
        if (c == -1)
            break;
        switch (c) {
        case 'd':
            options.debug = 1;
            break;
        case 'h':
            options.help = 1;
            break;
        case 'c':
            options.cutOff = atoi(optarg);
            break;
        case 'f':
            Feature = optarg;
            break;
        default:
            printf("?? getopt returned character code 0%o ??\n", c);
        }
    }
    /*
       if (optind == argc || options.help == 1) {
     */
    if (optind > argc || options.help == 1) {
        fprintf(stderr, "Usage: %s [options] [<] <SGA File>\n"
                "      - version %s\n"
                "      where options are:\n"
                "  \t\t -h     Show this help text\n"
                "  \t\t -d     Print debug information\n"
                "  \t\t -f     Feature specifications [name [+|-]]\n"
                "  \t\t -c     Count Cut-off (default is %d)\n"
                "\n\tThe program reads a ChIP-seq data file (or from stdin [<]) in SGA format (<SGA File>),\n"
                "\tand computes the total number of tag counts as well as the total sequence length.\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\n",
                argv[0], VERSION, options.cutOff);
        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);
    }
    /* Process Feature Specs */
    if (Feature == NULL) {
        ft_specs = 0;           /* Process all features */
    } else {
        options.ftName = 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;
            }
            options.ftName[i++] = *s++;
        }
        options.ftName[i] = '\0';
        while (isspace(*s++))
            options.strand = *s;
    }
    if (options.debug) {
        if (!ft_specs) {
            fprintf(stderr,
                    "Feature Specs: ALL -> Process all features\n");
        } else {
            if (options.strand == '\0')
                fprintf(stderr, "Feature Specs: Name : %s\n",
                        options.ftName);
            else
                fprintf(stderr, "Feature Specs: Name : %s  Strand %c\n",
                        options.ftName, options.strand);
        }
    }
    if (process_sga(input) != 0) {
        return 1;
    }
    if (input != stdin) {
        fclose(input);
    }
    return 0;
}
