/*
  filter_counts.c

  Filter Feature read counts that occur within user-defined regions
  of interest such as Repeat Mask regions.
  If the -k option is set, the program retains read counts that occur
  within user-defined regions of interest.

  # Arguments:
  # Feature name of the user-defined regions of interest (e.g. RMSK)


  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 retain;
} options_t;

static options_t options;

char *Feature = NULL;

int process_sga(FILE *input, char *iFile) {
    unsigned long pos;
    int cnt;
    char *s, *res, *buf;
    size_t bLen = LINE_SIZE;
    int prt = 0;
    char annotation[EXT_MAX];

    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 c = 1;
#endif
    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 = (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
        fprintf(stderr,
                " [%d] seq ID: %s   Feat: %s (%c)  Pos: %lu  Cnts: %d Ext: %s\n",
                c++, seq_id, ft, strand, pos, cnt, ext);
#endif
        if (!strcmp(ft, Feature)) {
            if (strand == '+') {
                prt++;
                strcpy(annotation, ext);
            } else if (strand == '-') {
                prt--;
            }
        } else {
            if (options.retain) {
                if (prt >= 1) {
                    if (!strcmp(ext, "")) {
                        printf("%s\t%s\t%lu\t%c\t%d\t%s\n", seq_id, ft,
                               pos, strand, cnt, annotation);
                    } else {
                        printf("%s\t%s\t%lu\t%c\t%d\t%s\t%s\n", seq_id, ft,
                               pos, strand, cnt, ext, annotation);
                    }
                }
            } else {
                if (prt < 1) {
                    if (!strcmp(ext, "")) {
                        printf("%s\t%s\t%lu\t%c\t%d\n", seq_id, ft, pos,
                               strand, cnt);
                    } else {
                        printf("%s\t%s\t%lu\t%c\t%d\t%s\n", seq_id, ft,
                               pos, strand, cnt, ext);
                    }
                }
            }
        }
    }                           /* End of While */
    free(s);
    return 0;
}

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

    while (1) {
        int c = getopt(argc, argv, "f:dhkr");
        if (c == -1)
            break;
        switch (c) {
        case 'f':
            Feature = optarg;
            break;
        case 'd':
            options.debug = 1;
            break;
        case 'h':
            options.help = 1;
            break;
        case 'r':
            options.retain = 1;
            break;
        case 'k':
            options.retain = 1;
            break;
        default:
            printf("?? getopt returned character code 0%o ??\n", c);
        }
    }
    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 defining the user-defined regions of interest (default=RMSK)\n"
                "  \t\t -k     Retain Mode on (\"keep\")\n"
                "  \t\t -r     Older synonym of -k (\"retain\")\n"
                "\n\tThe program reads a ChIP-seq data file (or from stdin [<]) in SGA format (<SGA File>)\n"
                "\tand filters out all read counts that occur within user-defined regions of interest, such\n"
                "\tas Repeat Mask regions (feature=RMSK).\n"
                "\tIf 'Retain Mode' is on, it only retains those lines that are within the user-defined regions,\n"
                "\tin which case the annotation of the user-defined regions is added to the output SGA file.\n\n",
                argv[0], VERSION);
        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, " Feature name: %s\n\n", Feature);
    }
    /* Process Feature */
    if (Feature == NULL) {
        if ((Feature = malloc(5 * sizeof(char))) == NULL) {
            perror("process_sga: malloc Feature");
            exit(1);
        }
        strcpy(Feature, "RMSK");
        free_flag = 1;
    } else {
        char *s = Feature;
        int i = 0;
        while (*s != 0 && !isspace(*s)) {
            if (i >= FT_MAX) {
                fprintf(stderr, "Feature Name too long \"%s\" \n",
                        Feature);
                return 1;
            }
            s++;
        }
    }
    if (process_sga(input, argv[optind++]) != 0) {
        return 1;
    }
    if (free_flag)
        free(Feature);
    return 0;
}
