/*==========================================================================
                SeqAn - The Library for Sequence Analysis
                          http://www.seqan.de 
  ============================================================================
  Copyright (C) 2010

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 3 of the License, or (at your option) any later version.

  This library 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
  Lesser General Public License for more details.
  ============================================================================
  Author: Manuel Holtgrewe <manuel.holtgrewe@fu-berlin.de>
  ============================================================================
  Read Analyzer code.

  Takes a reference sequence and Sam read alignment file as well as
  the sequencing technologies as its parameter and performs basic
  statistics on this file.  Currently, Illumina and 454 reads are
  supported.  For Illumina reads, the model used for the analysis is
  focused on edit distance whereas for 454, the focus is on explaining
  the difference in sequences by different homopolymer lengths.
  
  In order to work correctly, the Sam file put into this program
  should be generated by RazerS with the following options:

    -rr 100  # 100% recognition rate
    -i 90    # 10% error rate
    -of 4    # Sam output format
    -dr 0    # Only output the best matches for each read
    -id      # Allow indels
  ==========================================================================*/

#include <seqan/align.h>
#include <seqan/basic.h>
#include <seqan/store.h>

#include "read_analyzer.h"
#include "analyze_illumina.h"
#include "analyze_454.h"

using namespace seqan;


// Clear all but the contigs from the fragment store.
template <typename TFragmentStore>
void clearAllButContigs(TFragmentStore & fragmentStore)
{
  clear(fragmentStore.alignedReadStore);
  clear(fragmentStore.alignedReadTagStore);
  clear(fragmentStore.alignQualityStore);
  clear(fragmentStore.annotationNameStore);
  clear(fragmentStore.annotationStore);
  clear(fragmentStore.libraryNameStore);
  clear(fragmentStore.libraryStore);
  clear(fragmentStore.matePairNameStore);
  clear(fragmentStore.matePairStore);
  clear(fragmentStore.readNameStore);
  clear(fragmentStore.readSeqStore);
  clear(fragmentStore.readStore);
}


int main(int argc, char **argv) {
    // Check arguments.
    if (argc < 4 ||
        ((argc > 2) && (CharString(argv[1]) != "illumina" && CharString(argv[1]) != "454"))) {
        std::cerr << "Usage: read_analyzer illumina GENOME.FASTA FILE.Sam+" << std::endl
                  << "       read_analyzer 454 GENOME.FASTA FILE.Sam" << std::endl;
        return 1;
    }

    // =======================================================================
    // Load Contigs
    // =======================================================================
    typedef FragmentStore<> TFragmentStore;
    TFragmentStore fragmentStore;

    // Load contigs.
    std::cerr << "Reading FASTA contigs sequence file " << argv[2] << " ..." << std::endl;
    if (!loadContigs(fragmentStore, argv[2])) {
        std::cerr << "Could not read contigs." << std::endl;
        return 1;
    }

    // =======================================================================
    // Perform evaluation for each Sam file and print results.
    // =======================================================================
    if (CharString(argv[1]) == "illumina") {
        ReadEvaluationResult<Illumina> readResult;
        AlignmentEvaluationResult<Illumina> alignmentResult;
    
        for (int i = 3; i < argc; ++i) {
            clearAllButContigs(fragmentStore);

            // Read next Sam file.
            std::cerr << "Reading Sam file file " << argv[i] << " ..." << std::endl;
            std::fstream fstrm(argv[i], std::ios_base::in | std::ios_base::binary);
            if (!fstrm.is_open()) {
                std::cerr << "Could not open Sam file." << std::endl;
                return 1;
            }
            read(fstrm, fragmentStore, Sam());

            // Initialize counters.
            if (i == 3) {
              setReadLength(readResult, length(fragmentStore.readSeqStore[0]));
              setReadLength(alignmentResult, length(fragmentStore.readSeqStore[0]), length(fragmentStore.contigStore));
            }

            std::cerr << "Evaluating Illumina reads..." << std::endl;
            performReadEvaluation(readResult, fragmentStore);
            performAlignmentEvaluation(alignmentResult, fragmentStore);
        }

        // Print evaluation results.
        printReadEvaluationResults(readResult);
        printAlignmentEvaluationResults(alignmentResult, fragmentStore);
    } else {
        // Load Sam file.
        std::cerr << "Reading Sam file file " << argv[3] << " ..." << std::endl;
        std::fstream fstrm(argv[3], std::ios_base::in | std::ios_base::binary);
        if (!fstrm.is_open()) {
            std::cerr << "Could not open Sam file." << std::endl;
            return 1;
        }
        read(fstrm, fragmentStore, Sam());

        std::cerr << "Evaluating 454 reads..." << std::endl;

        ReadEvaluationResult<LS454> readResult;
        performReadEvaluation(readResult, fragmentStore);
        printReadEvaluationResults(readResult);
    }

    return 0;
}