Rcpp speed comparison for different ways of constructing sp objects

aaaline-tests.md

Compares Rcpp construction of sp::Line objects either:

1. Constructed as new objects with slots filled directly, or
2. Constructed the 'usual' way through filling the sp object with data.

In R terms, these are equivalent to:

library (sp)
# first way
obj <- new ("Line")
slot (obj, "coords") <- some_data
# second way
obj <- sp::Line (coords=some_data)

The actual Rcpp code constructs SpatialLines, with reasons for differences be readily seen from the sp source: direct construction of a SpatialLines object requires an sapply operation which has to be much slower than filling a new object, as the following code demonstrates.

First load the following C++ file and necessary R packages:

Sys.setenv ('PKG_CXXFLAGS'='-std=c++11')
Rcpp::sourceCpp('line-construction.cpp')
Sys.unsetenv ('PKG_CXXFLAGS')
library (microbenchmark)
library (sp)

Then the timing comparison:

mb_new <- microbenchmark ( obj <- line_new (), times=1000L )
tt_new <- formatC (mean (mb_new$time) / 1e6, format="f", digits=2) # in ms
mb_fill <- microbenchmark ( obj <- line_fill (), times=1000L )
tt_fill <- formatC (mean (mb_fill$time) / 1e6, format="f", digits=2)
cat ("construct new vs. direct fill times = ", tt_new, " vs ",
     tt_fill, "\n", sep="")

## construct new vs. direct fill times = 0.19 vs 0.96

line-construction.cpp

#include <string>
#include <fstream> // ifstream
#include <iostream>
#include <unordered_set>
#include <Rcpp.h>

using namespace Rcpp;

const float FLOAT_MAX = std::numeric_limits<float>::max ();

// [[Rcpp::export]]
S4 line_new ()
{
    const int nLines = 3;
    int nrow, count=0;
    float tempf, xmin = FLOAT_MAX, xmax = -FLOAT_MAX, 
          ymin = FLOAT_MAX, ymax = -FLOAT_MAX;
    std::vector <float> vec;
    List result (nLines);

    std::vector <std::string> colnames, rownames;
    colnames.push_back ("x");
    colnames.push_back ("y");
    List dimnames (0);

    Rcpp::Language line_call ("new", "Line");
    Rcpp::Language lines_call ("new", "Lines");
    Rcpp::S4 line;
    Rcpp::S4 lines;
    Rcpp::List dummy_list (0);

    for (int i=0; i<nLines; i++)
    {
        nrow = 10 - i * 2;
        rownames.resize (0);
        vec.resize (0);
        for (int j=0; j<(2 * nrow); j++)
            vec.push_back (10.0 * (float) i + (float) j);
        NumericMatrix nmat (Dimension (nrow, 2));
        for (int j=0; j<nrow; j++)
        {
            rownames.push_back ("r" + std::to_string (j));
            nmat (j, 0) = vec [j * 2];
            nmat (j, 1) = vec [j * 2 + 1];
            if (nmat (j, 0) < xmin)
                xmin = nmat (j, 0);
            else if (nmat (j, 0) > xmax)
                xmax = nmat (j, 0);
            if (nmat (j, 1) < ymin)
                ymin = nmat (j, 1);
            else if (nmat (j, 1) > ymax)
                ymax = nmat (j, 1);
        }
        dimnames.push_back (rownames);
        dimnames.push_back (colnames);
        nmat.attr ("dimnames") = dimnames;
        while (dimnames.size () > 0)
            dimnames.erase (0);

        line = line_call.eval ();
        line.slot ("coords") = nmat;
        dummy_list.push_back (line);
        lines = lines_call.eval ();
        lines.slot ("Lines") = dummy_list;
        lines.slot ("ID") = "number#" + std::to_string (i);
        result [count++] = lines;

        dummy_list.erase (0);
    }

    vec.resize (0);
    Rcpp::Language sp_lines_call ("new", "SpatialLines");
    Rcpp::S4 sp_lines;
    sp_lines = sp_lines_call.eval ();
    sp_lines.slot ("lines") = result;

    return sp_lines;
}

// [[Rcpp::export]]
S4 line_fill ()
{
    const int nLines = 3;
    int nrow, count=0;
    float tempf, xmin = FLOAT_MAX, xmax = -FLOAT_MAX, 
          ymin = FLOAT_MAX, ymax = -FLOAT_MAX;
    std::string tempstr;
    std::vector <float> vec;
    List result (nLines);

    std::vector <std::string> colnames, rownames;
    colnames.push_back ("x");
    colnames.push_back ("y");
    List dimnames (0);

    for (int i=0; i<nLines; i++)
    {
        nrow = 10 - i * 2;
        rownames.resize (0);
        vec.resize (0);
        for (int j=0; j<(2 * nrow); j++)
            vec.push_back (10.0 * (float) i + (float) j);
        NumericMatrix nmat (Dimension (nrow, 2));
        for (int j=0; j<nrow; j++)
        {
            rownames.push_back ("r" + std::to_string (j));
            nmat (j, 0) = vec [j * 2];
            nmat (j, 1) = vec [j * 2 + 1];
            if (nmat (j, 0) < xmin)
                xmin = nmat (j, 0);
            else if (nmat (j, 0) > xmax)
                xmax = nmat (j, 0);
            if (nmat (j, 1) < ymin)
                ymin = nmat (j, 1);
            else if (nmat (j, 1) > ymax)
                ymax = nmat (j, 1);
        }
        dimnames.push_back (rownames);
        dimnames.push_back (colnames);
        nmat.attr ("dimnames") = dimnames;
        while (dimnames.size () > 0)
            dimnames.erase (0);

        Rcpp::S4 line = Rcpp::Language ("Line", nmat).eval ();
        tempstr = "number#" + std::to_string (i);
        Rcpp::S4 lines = Rcpp::Language ("Lines", line, tempstr).eval ();
        result [count++] = lines;
    }

    vec.resize (0);
    Rcpp::Language sp_lines_call ("new", "SpatialLines");
    Rcpp::S4 sp_lines;
    sp_lines = sp_lines_call.eval ();
    sp_lines.slot ("lines") = result;

    return sp_lines;
}