Building a Simple Application Using 'libsequence'

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <stdexcept>
 
#include <Sequence/FST.hpp>
#include <Sequence/Fasta.hpp>
#include <Sequence/Alignment.hpp>
#include <Sequence/PolySites.hpp>
#include <Sequence/SeqExceptions.hpp>
 
/**
 * Converts from one simple streamable type to another.
 */
template <typename T, typename U>
T to_scalar(U from) {
    std::ostringstream o;
    o << from;
    std::istringstream i(o.str());
    T target;
    i >> target;
    if (i.fail() || !i.eof()) {
        throw std::runtime_error(o.str());
    }
    return target;
}
 
/**
 * Returns <code>true</code> if the second string is equal to the first n
 * characters of the first string, where n = length of the second string.
 *
 * @param s1        first string
 * @param s2        second string
 * @return          <code>true</code> if the first string begins with the
 *                  second
 */
bool startswith(const std::string& s1, const std::string& s2) {
    return s1.compare(0, s2.size(), s2) == 0;
}
 
/**
 * Parses command line invocation into a pair of vectors of strings,
 * representing vector of flag arguments (i.e., preceded by '-' or '--') and
 * positional arguments.
 */
std::pair<std::vector<std::string>, std::vector<std::string> > parse_cmd(
        int argc,
        char * argv[]) {
    std::vector<std::string> opts;
    std::vector<std::string> args;
    for (unsigned int i = 1; i < argc; ++i) {
        if (argv[i][0] == '-') {
            opts.push_back(argv[i]);
        } else {
            args.push_back(argv[i]);
        }
    }
    return std::make_pair(opts, args);
}
 
/**
 * Usage message.
 */
void show_usage(std::ostream& out=std::cerr) {
    out << "Usage: fst <FASTA-FILE> <POP-SIZE1> <POP-SIZE2> [POP-SIZE3 [POP-SIZE4 [...]]]\n";
}
 
/**
 * Usage message.
 */
void show_help(std::ostream& out=std::cout) {
    show_usage(out);
    out << std::endl;
    out << "Options:" << std::endl;
    out << "    --hsm1992      Use Hudson, Slatkin and Maddison (1992) formulation [default]" << std::endl;
    out << "    --slatkin1993  Use Slatkin (1993) formulation" << std::endl;
    out << "    --hbm1992      Use Hudson, Boos, and Kaplan (1992) formulation" << std::endl;
    out << "    --all          Use all of the above formulations (tab-delimited)" << std::endl;
}
 
/**
 * Main CLI and processing loop.
 */
int main(int argc, char * argv[]) {
 
    // Parsing and processing of invocation options and argument
 
    std::pair<std::vector<std::string>, std::vector<std::string> > cmd_opts = parse_cmd(argc, argv);
    std::vector<std::string> opts = cmd_opts.first;
    std::vector<std::string> args = cmd_opts.second;
    int method = 0;
 
    for (std::vector<std::string>::const_iterator oi = opts.begin(); oi != opts.end(); ++oi) {
        const std::string& opt = *oi;
        if (startswith(opt, "-h") || startswith(opt, "--he")) {
            show_help(std::cout);
            exit(0);
        } else if (startswith(opt, "--hsm1992")) {
            method = 0;
        } else if (startswith(opt, "--slatkin1993")) {
            method = 1;
        } else if (startswith(opt, "--hbm1992")) {
            method = 2;
        } else if (startswith(opt, "--all")) {
            method = -1;
        } else {
            std::cerr << "Unrecognized option '" << opt << "'. See '--help' for available options." << std::endl;
            exit(1);
        }
 
    }
 
    // Make sure mandatory arguments are specified
 
    if (args.size() < 3 ) {
        show_usage(std::cerr);
        exit(1);
    }
 
    // Input data filename
 
    std::string input_filepath = args[0];
 
    // Parsing and storing of subpopulation sizes
 
    std::vector<unsigned int> pop_sizes;
    pop_sizes.reserve(args.size() - 1);
    for (unsigned long i = 1; i < args.size(); ++i) {
        pop_sizes.push_back(to_scalar<unsigned int>(args[i]));
    }
 
    // Using libsequence to read data and convert to table of polymorphisms
 
    std::vector<Sequence::Fasta> alignment;
    Sequence::Alignment::GetData(alignment, input_filepath.c_str());
    Sequence::PolySites data(alignment);
 
    // Converting subpopulation sizes from vector to array
 
    unsigned int * config = &pop_sizes[0];
    Sequence::FST fst(&data, pop_sizes.size(), config);
 
    // Actual calls to libsequence's calculators and output
 
    if (method < 0) {
        std::cout << fst.HSM() << "\t" << fst.Slatkin() << "\t" << fst.HBK() << std::endl;
    } else if (method == 0) {
        std::cout << fst.HSM() << std::endl;
    } else if (method == 1) {
        std::cout << fst.Slatkin() << std::endl;
    } else if (method == 2) {
        std::cout << fst.HBK() << std::endl;
    }
}