5
\$\begingroup\$

I write quite a bit of code in C, but haven't done much C++ since my college CS classes. I have been revisiting C++ recently, and thought I would re-implement a program I had previously written in C, but this time in C++. I particularly took advantage of containers and other features standardized recently by the C++11 standard.

I really enjoyed writing the code, and have verified that the C++ version of my program generates the same output as the original C version. However, there is a substantial performance margin--I'm talking like an order of magnitude and then some for the test I ran.

Many programmers claim that C++ performance should be close to C. Of course there's a bit more overhead, but from 48 seconds to 18 minutes--that's just unreasonable.

I'm wondering whether I've done something completely wrong or inefficient, which is causing the performance discrepancy. Any feedback would be welcome.

Here is the C code:

/*

Copyright (c) 2013, Daniel S. Standage <[email protected]>

Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.

THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

*/

#include <getopt.h>
#include <stdio.h>
#include "khash.h"

//------------------------------------------------------------------------------
// Definitions/prototypes/initializations for data structures, functions, etc.
//------------------------------------------------------------------------------
#define MAX_LINE_LENGTH 8192
#define MAX_ID_LENGTH 1024
KHASH_MAP_INIT_STR(m32, unsigned)

typedef struct
{
  char delim;
  const char *outfile;
  FILE *outstream;
  unsigned numfiles;
} SmrOptions;

void smr_init_options(SmrOptions *options);
khash_t(m32) *smr_collect_molids(SmrOptions *options, khash_t(m32) **maps);
khash_t(m32) *smr_load_file(const char *filename);
void smr_parse_options(SmrOptions *options, int argc, char **argv);
void smr_print_matrix(SmrOptions *options, khash_t(m32) **maps);
void smr_print_usage(FILE *outstream);
void smr_terminate(SmrOptions *options, khash_t(m32) **maps);

//------------------------------------------------------------------------------
// Main method
//------------------------------------------------------------------------------

int main(int argc, char **argv)
{
  SmrOptions options;
  smr_init_options(&options);
  smr_parse_options(&options, argc, argv);

  unsigned i;
  khash_t(m32) **maps = malloc( sizeof(void *) * options.numfiles );
  for(i = 0; i < options.numfiles; i++)
    maps[i] = smr_load_file(argv[optind+i]);
  smr_print_matrix(&options, maps);

  smr_terminate(&options, maps);
  return 0;
}

//------------------------------------------------------------------------------
// Function implementations
//------------------------------------------------------------------------------
khash_t(m32) *smr_collect_molids(SmrOptions *options, khash_t(m32) **maps)
{
  unsigned i;
  khiter_t iter;
  khash_t(m32) *ids = kh_init(m32);
  for(i = 0; i < options->numfiles; i++)
  {
    for(iter = kh_begin(maps[i]); iter != kh_end(maps[i]); iter++)
    {
      if(!kh_exist(maps[i], iter))
        continue;
      const char *molid = kh_key(maps[i], iter);
      khint_t key = kh_get(m32, ids, molid);
      if(key == kh_end(ids))
      {
        int code;
        key = kh_put(m32, ids, molid, &code);
        if(!code)
        {
          fprintf(stderr, "error: failure storing key '%s'\n", molid);
          kh_del(m32, ids, key);
        }
        else
          kh_value(ids, key) = 1;
      }
    }
  }
  return ids;
}

void smr_init_options(SmrOptions *options)
{
  options->delim = ',';
  options->outfile = "stdout";
  options->outstream = stdout;
  options->numfiles = 0;
}

khash_t(m32) *smr_load_file(const char *filename)
{
  FILE *instream = fopen(filename, "r");
  if(instream == NULL)
  {
    fprintf(stderr, "error opening file '%s'\n", filename);
    exit(1);
  }

  khash_t(m32) *map = kh_init(m32);
  char buffer[MAX_LINE_LENGTH];
  while(fgets(buffer, MAX_LINE_LENGTH, instream) != NULL)
  {
    const char *tok = strtok(buffer, "\t\n");
    tok = strtok(NULL, "\t\n");
    int bflag = atoi(tok);
    if(bflag & 0x4)
      continue;

    tok = strtok(NULL, "\t\n");
    khint_t key = kh_get(m32, map, tok);

    if(key == kh_end(map))
    {
      int code;
      const char *molid = strdup(tok);
      key = kh_put(m32, map, molid, &code);
      if(!code)
      {
        fprintf(stderr, "error: failure storing key '%s'\n", molid);
        kh_del(m32, map, key);
      }
      else
        kh_value(map, key) = 0;
    }
    unsigned tsareadcount = kh_value(map, key);
    kh_value(map, key) = tsareadcount + 1;
  }

  fclose(instream);
  return map;
}

void smr_parse_options(SmrOptions *options, int argc, char **argv)
{
  int opt = 0;
  int optindex = 0;
  const char *optstr = "d:ho:";
  const struct option smr_options[] =
  {
    { "delim", required_argument, NULL, 'd' },
    { "help", no_argument, NULL, 'h' },
    { "outfile", required_argument, NULL, 'o' },
    { NULL, no_argument, NULL, 0 },
  };

  for(opt = getopt_long(argc, argv, optstr, smr_options, &optindex);
      opt != -1;
      opt = getopt_long(argc, argv, optstr, smr_options, &optindex))
  {
    switch(opt)
    {
      case 'd':
        if(strcmp(optarg, "\\t") == 0)
          optarg = "\t";
        else if(strlen(optarg) > 1)
        {
          fprintf(stderr, "warning: string '%s' provided for delimiter, using "
                  "only '%c'\n", optarg, optarg[0]);
        }
        options->delim = optarg[0];
        break;
      case 'h':
        smr_print_usage(stdout);
        exit(0);
        break;
      case 'o':
        options->outfile = optarg;
        break;
      default:
        fprintf(stderr, "error: unknown option '%c'\n", opt);
        smr_print_usage(stderr);
        break;
    }
  }

  if(strcmp(options->outfile, "stdout") != 0)
  {
    options->outstream = fopen(options->outfile, "w");
    if(options->outstream == NULL)
    {
      fprintf(stderr, "error: unable to open output file '%s'\n",
              options->outfile);
      exit(1);
    }
  }

  options->numfiles = argc - optind;
  if(options->numfiles < 1)
  {
    fputs("expected 1 or more input files\n", stderr);
    smr_print_usage(stderr);
    exit(1);
  }
}

void smr_print_matrix(SmrOptions *options, khash_t(m32) **maps)
{
  khiter_t iter;
  unsigned i;
  khash_t(m32) *molids = smr_collect_molids(options, maps);
  for(iter = kh_begin(molids); iter != kh_end(molids); iter++)
  {
    if(!kh_exist(molids, iter))
      continue;
    const char *molid = kh_key(molids, iter);

    fprintf(options->outstream, "%s%c", molid, options->delim);
    for(i = 0; i < options->numfiles; i++)
    {
      if(i > 0)
        fputc(options->delim, options->outstream);

      khint_t key = kh_get(m32, maps[i], molid);
      if(key == kh_end(maps[i]))
        fputc('0', options->outstream);
      else
      {
        unsigned readcount = kh_value(maps[i], key);
        fprintf(options->outstream, "%u", readcount);
      }
    }
    fputc('\n', options->outstream);
  }
  kh_destroy(m32, molids);
}

void smr_print_usage(FILE *outstream)
{
  fputs("\nSMR: SAM mapped reads\n\n"
"The input to SMR is 1 or more SAM files. The output is a table (1 column for\n"
"each input file) showing the number of reads that map to each sequence.\n\n"
"Usage: smr [options] sample-1.sam sample-2.sam ... sample-n.sam\n"
" Options:\n"
" -d|--delim: CHAR delimiter for output data; default is comma\n"
" -h|--help print this help message and exit\n"
" -o|--outfile: FILE name of file to which read counts will be\n"
" written; default is terminal (stdout)\n",
        outstream);
}

void smr_terminate(SmrOptions *options, khash_t(m32) **maps)
{
  unsigned i;
  khint_t iter;
  for(i = 0; i < options->numfiles; i++)
  {
    for(iter = kh_begin(maps[i]); iter != kh_end(maps[i]); iter++)
    {
      if(!kh_exist(maps[i], iter))
        continue;
      char *molid = (char *)kh_key(maps[i], iter);
      free(molid);
    }
    kh_destroy(m32, maps[i]);
  }
  free(maps);
  fclose(options->outstream);
}

And here is the C++ code:

/*

Copyright (c) 2013, Daniel S. Standage <[email protected]>

Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.

THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

*/

#include <unistd.h>
#include <cstring>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>

typedef std::unordered_set<std::string> uset;
typedef std::unordered_map<std::string, unsigned> umap;
struct SmrOptions
{
  char delim;
  std::string outfile;
  std::ostream& outstream;
  unsigned numfiles;
  std::vector<std::string> infiles;
  SmrOptions(char d, std::string of, std::ostream& os, unsigned nf,
             std::vector<std::string> infls)
  : delim(d), outfile(of), outstream(os), numfiles(nf), infiles(infls) {}
};
typedef struct SmrOptions SmrOptions;

uset smr_collect_molids(SmrOptions& options, std::vector<umap>& rm2seqPerSample);
umap smr_load_file(std::istream& instream, char delim);
SmrOptions smr_parse_options(int argc, char **argv);
void smr_print_matrix(SmrOptions& options, std::vector<umap>& rm2seqPerSample);
void smr_print_usage(std::ostream& outstream);
std::vector<std::string>& smr_string_split(const std::string &s, char delim,
                                           std::vector<std::string> &elems);

//------------------------------------------------------------------------------
// Main method
//------------------------------------------------------------------------------

int main(int argc, char **argv)
{
  SmrOptions options = smr_parse_options(argc, argv);

  std::vector<umap> rm2seqPerSample;
  for(unsigned i = 0; i < options.numfiles; i++)
  {
    std::ifstream ifs (std::string(options.infiles[i]));
    if(!ifs.is_open())
    {
      std::cerr << "error opening input file "
                << options.infiles[i]
                << std::endl;
      exit(1);
    }
    rm2seqPerSample.push_back(smr_load_file(ifs, options.delim));
    ifs.close();
  }
  smr_print_matrix(options, rm2seqPerSample);

  //options.outstream.close();
  return 0;
}

//------------------------------------------------------------------------------
// Function implementations
//------------------------------------------------------------------------------
uset smr_collect_molids(SmrOptions& options, std::vector<umap>& rm2seqPerSample)
{
  uset molids;
  for(std::vector<umap>::iterator ssiter = rm2seqPerSample.begin();
      ssiter != rm2seqPerSample.end();
      ssiter++)
  {
    umap& rm2seq = *ssiter;
    for(umap::iterator iter = rm2seq.begin(); iter != rm2seq.end(); iter++)
    {
      molids.emplace(iter->first);
    }
  }
  return molids;
}

umap smr_load_file(std::istream& instream, char delim)
{
  umap rm2seq;

  std::string buffer;
  while(std::getline(instream, buffer))
  {
    if(buffer[0] == '@')
      continue;

    std::vector<std::string> tokens;
    smr_string_split(buffer, '\t', tokens);
    std::string molid = tokens[2];
    std::string bflag_str = tokens[1];
    int bflag = std::stoi(bflag_str);
    if(bflag & 0x4)
      continue;

    umap::iterator keyvaluepair = rm2seq.find(molid);
    if(keyvaluepair == rm2seq.end())
      rm2seq.emplace(molid, 1);
    else
      rm2seq[molid] += 1;
  }

  return rm2seq;
}

SmrOptions smr_parse_options(int argc, char **argv)
{
  char delim = ',';
  std::string outfile = "stdout";

  char opt;
  const char *arg;
  while((opt = getopt(argc, argv, "d:ho:")) != -1)
  {
    switch(opt)
    {
      case 'd':
        arg = optarg;
        if(strcmp(optarg, "\\t") == 0)
          arg = "\t";
        else if(strlen(optarg) > 1)
        {
          std::cerr << "warning: string '"
                    << arg
                    << "' provided for delimiter, using only '"
                    << arg[0]
                    << "'"
                    << std::endl;
        }
        delim = arg[0];
        break;
      case 'h':
        smr_print_usage(std::cout);
        exit(0);
        break;
      case 'o':
        outfile = optarg;
        break;
      default:
        fprintf(stderr, "error: unknown option '%c'\n", opt);
        smr_print_usage(std::cerr);
        break;
    }
  }

  std::ofstream outfilestream;
  bool usestdout = (outfile.compare(std::string("stdout")) == 0);
  if(!usestdout)
    outfilestream.open(outfile, std::ios::out);
  std::ostream& outstream = (usestdout ? std::cout : outfilestream);

  unsigned numfiles = argc - optind;
  if(numfiles < 1)
  {
    std::cerr << "expected 1 or more input files" << std::endl;
    smr_print_usage(std::cerr);
    exit(1);
  }

  std::vector<std::string> infiles;
  for(unsigned i = 0; i < numfiles; i++)
  {
    int ind = optind + i;
    std::string filename (argv[ind]);
    infiles.push_back(filename);
  }

  return SmrOptions (delim, outfile, outstream, numfiles, infiles);
}

void smr_print_matrix(SmrOptions& options, std::vector<umap>& rm2seqPerSample)
{
  uset molids = smr_collect_molids(options, rm2seqPerSample);
  for(uset::iterator iter = molids.begin(); iter != molids.end(); iter++)
  {
    std::string molid = *iter;
    options.outstream << molid << options.delim;

    for(std::vector<umap>::iterator sampleiter = rm2seqPerSample.begin();
        sampleiter != rm2seqPerSample.end();
        sampleiter++)
    {
      if(sampleiter != rm2seqPerSample.begin())
        options.outstream << options.delim;

      umap rm2seq = *sampleiter;
      umap::const_iterator keyvaluepair = rm2seq.find(molid);
      if(keyvaluepair == rm2seq.end())
        options.outstream << '0';
      else
        options.outstream << keyvaluepair->second;
    }
    options.outstream << std::endl;
  }
}

void smr_print_usage(std::ostream& outstream)
{
  outstream << std::endl << "SMR: SAM mapped reads" << std::endl << std::endl
            << "The input to SMR is 1 or more SAM files. The output is a table (1 column for" << std::endl
            << "each input file) showing the number of reads that map to each sequence." << std::endl << std::endl
            << "Usage: smr [options] sample-1.sam sample-2.sam ... sample-n.sam" << std::endl
            << " Options:" << std::endl
            << " -d|--delim: CHAR delimiter for output data; default is comma" << std::endl
            << " -h|--help print this help message and exit" << std::endl
            << " -o|--outfile: FILE name of file to which read counts will be" << std::endl
            << " written; default is terminal (stdout)" << std::endl << std::endl;
}

std::vector<std::string>& smr_string_split(const std::string &s, char delim,
                                           std::vector<std::string> &elems)
{
    std::stringstream ss(s);
    std::string item;
    while(std::getline(ss, item, delim))
    {
      elems.push_back(item);
    }
    return elems;
}

These are accessible from commit c941a4aee4 of my GitHub repo.

\$\endgroup\$
10
  • 2
    \$\begingroup\$ Lot of extra copying. But it's hard to measure actual throughput with a set of input data. Also to help other people in the future we ask you post the code for review on the site. This has two effects: 1) It makes sure the code that we spend time reviewing is available for others to compare review notes against 2) The code is in the public domain (see last line on page). \$\endgroup\$ Nov 2, 2013 at 7:16
  • \$\begingroup\$ @LokiAstari Did you mean hard to measure without input data? A sample input file is available at de.iplantcollaborative.org/dl/…, but be warned its 26GB. I'm working on uploading a compressed file, will update once I have. \$\endgroup\$ Nov 2, 2013 at 11:33
  • \$\begingroup\$ Ok, here is the link to the same file gzip-compressed: de.iplantcollaborative.org/dl/…. Any subset of this file will work fine for testing each line is an independent data entry. \$\endgroup\$ Nov 2, 2013 at 12:36
  • 2
    \$\begingroup\$ @LokiAstari maybe nitpicking, but creative commons != public domain. We have some restrictions such as attribution required and share alike. \$\endgroup\$
    – Adam
    Nov 2, 2013 at 23:47
  • \$\begingroup\$ Who said C is faster than C++ at all? stackoverflow.com/questions/6955114/is-c-notably-faster-than-c \$\endgroup\$
    – Celeritas
    Nov 3, 2013 at 2:26

1 Answer 1

5
+50
\$\begingroup\$

As the comments alluded to, there's a lot of copying going on. Since I don't have profile data (wrong way to go after performance), I'll go over differences in the approach.

main()

In main(), the C version pre-allocates an array of options.numfiles pointers, but the C++ version merely uses push_back() to build a vector. Since you know the size, you could make the C++ version pre-allocate similarly by calling rm2seqPerSample.reserve(options.numfiles).

smr_load_file returns a pointer in the C version, and a std::unordered_map<string, unsigned> in the C++ version. While this probably uses RVO to avoid an initial copy or move, this is then pushed into the above vector, which must move or copy the map. It would be more idiomatic in C++ to wrap this into a class, and perhaps call emplace_back instead:

struct smr_file : private std::unordered_map<std::string, unsigned>
{
    smr_file(std::istream& instream, char delim)
    {
        // code from smr_load_file operating on *this here
    }
};
: : :
int main(...)
{
    std::vector<smr_file> rm2secPerSample;
    : : :
    rm2seqPerSample.emplace_back(ifs, options.delim);
    : : :
}

smr_load_file()

Diving into the smr_load_file, there are small things and larger things. Starting small, you can replace rm2seq[molid] += 1 with keyvaluepair->second += 1 to avoid a second lookup. Or maybe you can replace the entire find/emplace/increment stanza with just rm2seq[molid] += 1, if you can ensure the key-value-pair instantiates with value 0. (The former could help with efficiencies; the latter could help with local code complexity, which is usually more important.) Going a bit larger, the call to smr_string_split worried me, as string splits are easy to get wrong, and it sounds like this could well be one of your bottlenecks.

smr_string_split()

Looking at smr_string_split, there are a few things I would likely change. First off, a badly formatted line could result in a crash in the calling code (if there's fewer than two delimiters found, only one or two elements will be in the vector, so the following tokens[2] will misbehave). Second, there are a lot of data copies in these few lines of code, such as into the stringstream, out into the string, and onward into the vector<string> elems. I would be tempted to replace this with a more straightforward method that specifically parses out these items, either returning through "out" parameters (bool smr_parse_line(const std::string& line, char delim, int &bflag, std::string& molid)), with a tuple (std::tuple<int, std::string> smr_parse_line(const std::string& line, char delim)), or even just implement it inline. As for implementing this method, I would tend towards std::string::find and substr calls, as it seems you only support a predefined number of tokens.

smr_print_matrix()

smr_print_matrix starts with a copy or move, though probably elided due to RVO, but then loops suboptimally. The manual iterator bounds checking in the iterator faces a potential performance penalty due to its repeated calculation of end(). I would suggest replacing the for loops with a range for: for (auto& molid : molids), similarly for (auto& rm2seq : rm2seqPerSample). The only trick is handling your conditional output of options.delim. Here your use of rm2seq.find is correct, although I can't guess whether conditionally outputting a char instead of an int would have any effect. If it doesn't, something like options.outstream << (keyvaluepair != rm2seq.end()) ? keyvaluepair->second : 0; might be a maintenance win (or loss if you hate ?:).

smr_collect_molids()

The for loops should be changed similarly, but at least here you use references to capture your iterator's referenced value. The translation between C and C++ here seems a little iffy, as I don't see what captures the else kh_value(ids, key) = 1.

overview

Your code has taken a first stab at becoming idiomatic C++ code, but it has more steps to go. You've replaced the khash data structure with data structures from C++, but are still using them in a largely C fashion. There's nothing inherently wrong with this, although it may surprise many C++ developers, and may be unable to leverage many of the improvements in C++.

The biggest general thing to watch out for are invisible copies of large amounts of data. Some of this can be dodged by making certain to call methods that support move. Other parts can be removed by changing the ownership of objects to only create them in their final place (such as my smr_file example). Still other parts can be removed by not changing between different representations of the same data (string to stringstream, and so forth). Finally, just because it is C++, don't feel you have to give up C tools that happen to be faster at certain tasks. Instead try to understand why they are faster, and see if there are alternatives in C++ that you didn't already consider.

\$\endgroup\$
3
  • \$\begingroup\$ Thanks for your feedback. The pointers regarding idiomatic C++ were particularly helpful. Most of this has not made any difference in performance, unfortunately. It wasn't until I addressed the smr_string_split function (gave up on stringstream and reverted to strtok) that I was able to achieve C-like performance. \$\endgroup\$ Nov 14, 2013 at 3:22
  • \$\begingroup\$ I'm not surprised that smr_string_split was the bottleneck. I don't want to fan the iostream flames, but all the copies there just couldn't help. By the way, in theory string::find and string::substr should be faster than strtok, as they throw away and recalculate less information (length of substring), but it's hard to guess if the difference would be measurable. \$\endgroup\$ Nov 14, 2013 at 3:32
  • \$\begingroup\$ Actually, now that I've refactored a bit more, it seems that the string parsing was only a small bit of the performance issue. The biggest performance issue seems to have been the use of the C++ input and output streams. Replacing these with C-style file pointers (which I had done to some extent while fixing the string parsing code) is what really makes the biggest difference. \$\endgroup\$ Nov 14, 2013 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.