/*
Utility to extract a given subset of loci from each of a set of pileup files,
returning the sorted pileup lines.

Compile with:
g++ -O3 -DNDEBUG pileup_loci_comparison.cc -o /full/path/to/pileup_loci_comparison

Use:
/full/path/to/pileup_loci_comparison
*/

#include <cstdio>
#include <libgen.h>

#include <vector>
#include <set>
#include <string>


int main(int argc, char ** argv)
{
    FILE * pileup_fh = fopen(argv[1], "r");
    FILE * loci_fh = fopen(argv[2], "r");
    FILE * output_fh = fopen(argv[3],"w");

    if (argc != 4)
    {
        printf("\nusage: pileup_loci_comparison <pileup_file_list> <loci_list> <output_file>\n"
               "\n\n"
               "<pileup_file_list> is a file with full pathnames to each pileup file, one per line\n\n"
               "<loci_list> is a file with one locus per line formatted as <chromosome>[tab]<position>\n\n"
               "<output_file> will contain the selected pileup lines in order provided by <pileup_file_list> and <loci_list>\n\n"
               "\n"
               "Note: <loci_list> must be sorted in chromosome order and locus order matching all pileup files\n\n");
        exit(0);
    }

    if (pileup_fh == NULL)
    {
        fprintf(stderr, "Couldn't open pileup list file %s\n", argv[1]);
        exit(1);
    }
    if (loci_fh == NULL)
    {
        fprintf(stderr, "Couldn't open loci file %s\n", argv[2]);
        exit(1);
    }
    if (output_fh == NULL)
    {
        fprintf(stderr, "Couldn't open output file %s\n", argv[3]);
    }

    char chromosome[1024];
    int locus;

    std::set<std::pair<std::string, int> > loci;


    while (! feof(loci_fh))
    {
        fscanf(loci_fh, "%s\t%i\n", chromosome, &locus);
        loci.insert(std::make_pair(std::string(chromosome), locus));
    }
    fclose(loci_fh);

    char pileup_file[256];

    FILE * pileup_file_fh;
    std::vector<FILE *> pileup_fhs;
    std::vector<std::string> pileup_files;

    while (! feof(pileup_fh))
    {
        fscanf(pileup_fh, "%s\n", pileup_file);
        pileup_file_fh = fopen(pileup_file, "r");
        if (pileup_file_fh == NULL)
        {
            fprintf(stderr, "Couldn't open pileup file %s\n", pileup_file);
            exit(1);
        }

        pileup_fhs.push_back(pileup_file_fh);
        pileup_files.push_back(std::string(basename(pileup_file)));
    }
    fclose(pileup_fh);

    char mutant[1024];
    char refbase;
    int depth;
    char * pileup = new char[1000000];
    std::pair<std::string, int> next_locus;

    while (! loci.empty())
    {
        next_locus = *loci.begin();
        loci.erase(loci.begin());

        for (size_t p = 0; p != pileup_fhs.size(); ++p)
        {

            while(1)
            {
                fscanf(pileup_fhs[p], "%s\t%i\t%c\t%i\t%s\n",
                       chromosome, &locus, &refbase, &depth, pileup);
              
                if (strcmp(chromosome, next_locus.first.c_str()) == 0
                    && locus >= next_locus.second)
                {
                    //we are at or past the locus of interest
                    if (locus == next_locus.second)
                    {
                        fprintf(output_fh, "%s\t%s\t%i\t%c\t%i\t%s\n",
                                pileup_files[p].c_str(), chromosome, locus, refbase, depth, pileup);
                    }
                    break;
                }
            }
        }
        fflush(output_fh);
    }

    for (size_t p = 0; p != pileup_fhs.size(); ++p)
    {
        fclose(pileup_fhs[p]);
    }

    fclose(output_fh);
    delete pileup;

    return 0;
}