Using RcppProgress to control the long computations in C++

Karl Forner — written May 16, 2013 — source

Usually you write C++ code with R when you want to speedup some calculations. Depending on the parameters, and especially during the development, it is difficult to anticipate the execution time of your computation, so that you do not know if you have to wait for one minute or several hours.

RcppProgress is a tool to help you monitor the execution time of your C++ code, by providing a way to interrupt the execution inside the C++ code, and also to display a progress bar indicative of the state of your computation.

Additionally, it is compatible with multithreaded code, for example using OpenMP, which is not as trivial as it may seem since you cannot just stop the execution in one thread. Also, not all threads should be writing in the console to avoid garbled output.

// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
// [[Rcpp::export]]
double long_computation(int nb) {
    double sum = 0;
    for (int i = 0; i < nb; ++i) {
        for (int j = 0; j < nb; ++j) {
	    sum += R::dlnorm(i+j, 0.0, 1.0, 0);
	}
    }
    return sum + nb;
}
    system.time(s  <- long_computation(1000))
   user  system elapsed 
  0.116   0.000   0.116 
    s
[1] 1002

Checking for user interrupts

Let’s modify our code to add a check for user interruption by calling the function Progress::check_abort. Note the Rcpp::depends(RcppProgress) attribute in the header part that takes care of the include path for the progress.hpp header.

Now the long_computation2 call should be interruptible (with CTRL+C in the classic R console).

// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
// [[Rcpp::export]]
double long_computation2(int nb) {
    double sum = 0;
    Progress p(0, false); // we need an instance, should be improved in next version
    for (int i = 0; i < nb; ++i) {
        if (Progress::check_abort() )
            return -1.0;
  	for (int j = 0; j < nb; ++j) {
	    sum += R::dlnorm(i+j, 0.0, 1.0, 0);
	}
    }
    return sum + nb;
}
    system.time(s  <- long_computation2(3000)) # interrupt me
   user  system elapsed 
  1.044   0.000   1.044 
    s
[1] 3002

You may wonder why we put the check_abort call in the first loop instead that in the second. The performance cost of check_abort call is not negligible. It should be put in a place called often enough (once per second) yet not too often.

Adding a progress bar

Time to add the progress bar. The increment function is quite fast, so we can put it in the second loop. In real life example, it is sufficient to put it at a place called at least every second.

// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
// [[Rcpp::export]]
double long_computation3(int nb, bool display_progress=true) {
    double sum = 0;
    Progress p(nb*nb, display_progress);
    for (int i = 0; i < nb; ++i) {
        if (Progress::check_abort() )
            return -1.0;
        for (int j = 0; j < nb; ++j) {
            p.increment(); // update progress
	    sum += R::dlnorm(i+j, 0.0, 1.0, 0);
	}
    }
    return sum + nb;
}
    system.time(s  <- long_computation3(3000)) # interrupt me
   user  system elapsed 
  1.072   0.000   1.073 
    s
[1] 3002

OpenMP support

First we need this to enable OpenMP support for gcc:

Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
Sys.setenv("PKG_LIBS"="-fopenmp")

Future Rcpp versions should have a plugin which does this for us.

Here is an OpenMP version of our function:

#ifdef _OPENMP
#include <omp.h>
#endif
// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
// [[Rcpp::export]]
double long_computation_omp(int nb, int threads=1) {
#ifdef _OPENMP
    if ( threads > 0 )
        omp_set_num_threads( threads );
    REprintf("Number of threads=%i\n", omp_get_max_threads());
#endif
 
    double sum = 0;
#pragma omp parallel for schedule(dynamic)   
    for (int i = 0; i < nb; ++i) {
        double thread_sum = 0;
  	for (int j = 0; j < nb; ++j) {
	    thread_sum += R::dlnorm(i+j, 0.0, 1.0, 0);
	}
        sum += thread_sum;
    }
    return sum + nb;
}

Now check that it is parallelized:

    system.time(s4 <- long_computation_omp(5000, 4))
   user  system elapsed 
  4.208   0.000   1.069 
    s4
[1] 5002
    system.time(s1 <- long_computation_omp(5000, 1))
   user  system elapsed 
  2.948   0.000   2.949 
    s1
[1] 5002

adding progress monitoring to the openMP function

#ifdef _OPENMP
#include <omp.h>
#endif
// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
// [[Rcpp::export]]
double long_computation_omp2(int nb, int threads=1) {
#ifdef _OPENMP
    if ( threads > 0 )
        omp_set_num_threads( threads );
#endif
    Progress p(nb, true);
    double sum = 0;
#pragma omp parallel for schedule(dynamic)   
    for (int i = 0; i < nb; ++i) {
        double thread_sum = 0;
        if ( ! Progress::check_abort() ) {
            p.increment(); // update progress
            for (int j = 0; j < nb; ++j) {
                thread_sum += R::dlnorm(i+j, 0.0, 1.0, 0);
            }
        }
        sum += thread_sum;
    }
    return sum + nb;
}
    system.time(s <- long_computation_omp2(5000, 4))
   user  system elapsed 
  4.108   0.000   1.049 

Test it now

If you want to test it now in your R console, just paste the following code (after installing the RcppProgress package, of course):

library(Rcpp)
Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
Sys.setenv("PKG_LIBS"="-fopenmp")

code='
#ifdef _OPENMP
#include <omp.h>
#endif
// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>

// [[Rcpp::export]]
double long_computation_omp2(int nb, int threads=1) {
#ifdef _OPENMP
    if ( threads > 0 )
        omp_set_num_threads( threads );
    REprintf("Number of threads=%i\\n", omp_get_max_threads());
#endif
    Progress p(nb, true);
    double sum = 0;
#pragma omp parallel for schedule(dynamic)   
    for (int i = 0; i < nb; ++i) {
        double thread_sum = 0;
        if ( ! Progress::check_abort() ) {
            p.increment(); // update progress
            for (int j = 0; j < nb; ++j) {
                thread_sum += R::dlnorm(i+j, 0.0, 1.0, 0);
            }
        }
        sum += thread_sum;
    }
  
    return sum + nb;
}
'

sourceCpp(code=code)
s <- long_computation_omp2(10000, 4)

Karl Forner
Quartz Bio

tags: interrupt  openmp 

Related Articles