Rcpp: reading file functions codes

/*
All these function are useful in Reading file.read_file_delime function read the file from start to end.'Tol' is an option to make string into lower case.
                              MACHERKI M E.
                              16/02/2016


Download source file:link

*/


#include <fstream>
#include <sstream>
#include <string>
#include <algorithm>
#include<cctype>
#include<Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
CharacterVector read_file(std::string path){
std::ifstream t(path.c_str());//connect with file 
std::stringstream ss;
ss<<t.rdbuf();                // scan file
return ss.str();
}
//[[Rcpp::export]]
CharacterVector read_file1(std::string path){
std::ifstream in(path.c_str());
std::string contents;
in.seekg(0,std::ios::end);
contents.resize(in.tellg());
in.seekg(0,std::ios::beg);
in.read(&contents[0],contents.size());
in.close();
return contents;
}
// function to read file from start to end pointer
//[[Rcpp::export]]
CharacterVector read_file_delime(std::string path,int start,int end,int TOL){
std::ifstream in(path.c_str());
std::string contents;
in.seekg(0,std::ios::end);
contents.resize(in.tellg());
in.seekg(start,std::ios::beg);
in.read(&contents[0],end);
in.close();
contents.erase(std::remove(contents.begin(),contents.end(),'\n'),contents.end());
//remove lines delimiter
if(TOL>0){//force to lower case transformation
std::transform(contents.begin(),contents.end(),contents.begin(),tolower);
}
return contents.substr(start-1,end-start);
}
/* read a file from start  to end*/
/* set the file path before using */

/*** R
###*******************choose file to read************************####
con<-utils::choose.files()
system.time(FR<-readLines(con))     #vector as the split(x,'\n') 
system.time(Fcpp1<-read_file(con))  #contents Line delimiter
system.time(Fcpp2<-read_file1(con)) #contents Line delimiter
#one string without \n and in lower case
system.time(Fcpp3<-read_file_delime(con,20,100,1)) 

*/ 

Merge sort: Rcpp integration code

/*sorting function using the merge sort algorithm
*/

#include <Rcpp.h>
using namespace Rcpp;

void intercal(int p, int q, int r, NumericVector v,  NumericVector w)
{
  int i, j, k;
  i = p;
  j = q;
  k = 0;
  while (i < q && j < r) {
    if (v[i] < v[j]) {
      w[k] = v[i];
      i++;
    }
    else {
      w[k] = v[j];
      j++;
    }
    k++;
  }
  while (i < q) {
    w[k] = v[i];
    i++;
    k++;
  }
  while (j < r) {
    w[k] = v[j];
    j++;
    k++;
  }
  for (i = p; i < r; i++)
    v[i] = w[i-p];
}

void mergesort(int p, int r, NumericVector v, NumericVector aux)
{
  int q;
  if (p < r - 1) {
    q = (p + r) / 2;
    mergesort(p, q, v,aux);
    mergesort(q, r, v,aux);
    intercal(p, q, r, v,aux);
  }
}

// [[Rcpp::export]]
NumericVector sort_merge(NumericVector vetor) {
  Rcpp::NumericVector res = Rcpp::clone(vetor);
  Rcpp::NumericVector aux = Rcpp::clone(vetor);
  int n = res.size();
  mergesort(0,n,res,aux);
  return res;}

The first code:Rcpp

sourceCpp is the important function in Rcpp. The Rcpp code is written into file with'.cpp' extension. The cppFunction can be useful but we can not  easily combine many functions in the same time. This is an example to use the two possibility:

• ccpFunction

R code:

Rcpp::cppFunction(
'int fibonacci(const int x) {
if (x == 0) return(0);
if (x == 1) return(1);
return (fibonacci(x - 1)) + fibonacci(x - 2);
}')

fibonacci(20)

• sourceCpp

I used the example used in the Rstudio C++ source file.

Get the source file: firstrun.cpp and save it in the same R directory.

R code:

Rcpp::sourceCpp('firstrun.cpp')

Rcpp:the first step

Rcpp uses the C++ compiler to produce a sharing library connectable with R environment. The fist step is to install  the Rtools program provided in CRAN site. The second step is to install the Rcpp package.It is important to use a good  editor such Rstudio which integrates the C++ code directly.

Why I use Rcpp ?

Rcpp is an R add-on package which facilitates extending R with C++ functions. It is being used for anything from small and quickly constructed add-on functions written either to fluidly experiment with something new or to accelerate computing by replacing an R function with its C++ equivalent to large-scale bindings for existing libraries, or as a building block in entirely new research computing environments. While still relatively new as a project, Rcpp has already become widely deployed among users and developers in the R community. Rcpp is now the most popular language extension for the R system and used by over 500 CRAN packages as well as ten BioConductor packages. This site aims to provide a solid introduction to Rcpp and provides many codes of the most common algorithms.

Numerical coding of DNA: Rcpp code

/*
In general, DNA sequences are  stored as character data into different formats
of files like FASTA,FNA... ALso numeric coding of DNA in base 4 exists  such 
'A'=0,'G'=1,'C'=2 and 'T'=3. But why I use this kind of coding? 
The first step of sequence alignment algorithms such BLAST or FASTA consists on 
stepwise splitting of the DNA sequence into oligos of 11 bases. This operation 
influences  the execution  running time of the program.
This is an example of codes  ;the first use string manipulation and the second 
use numeric coding in base 5(to overcame  N in the sequence).
PDF file link:source file


                     MACHERKI M.E
                     14/02/2016
 */
 
 /*
 This source file  uses C++ code.You have to install Rcpp package to run it.
 A FASTA file  path is required to test this code in R.Please check the bottom 
 of the page.
 */ 
 
#include <fstream>
#include<Rcpp.h>
#include<string.h>
#include <sstream>
#include<cctype>
#include <algorithm>
using namespace Rcpp;

/*
First of all, this is a function to read FASTA file into R
-----!!!!!!!!!!!!! one sequence only!!!!!!!!!!!!!!-----------
*/
//[[Rcpp::export]]
inline CharacterVector read_fasta(std::string file ){
CharacterVector records;
std::ifstream in(file.c_str());
in.get();// remove first >
std::string rec ;
while(getline(in,rec,'>')){
int newLineLoc=rec.find('\n');
std::string header =rec.substr(0,newLineLoc);
std::string sequence= rec.substr(newLineLoc+1,rec.length()-newLineLoc-2);
sequence.erase(std::remove(sequence.begin(),sequence.end(),'\n'),sequence.end());
std::transform(sequence.begin(),sequence.end(),sequence.begin(),tolower);
records[header]=sequence;
}
return records;
}
/*
String splitter function into N size oligos (Only for lower case string).

*/
// [[Rcpp::export]]
 inline CharacterVector SPLIT(const std::string  strings,const int n){
    int lim=strings.length()-n+1 ;//size of sub string vector
    CharacterVector resultat(lim) ;
          for ( int j=0; j<lim ;j++) resultat[j] = strings.substr(j,n) ;
return resultat ;
} 
/*
String splitter function into N size oligos using numeric coding.
I used a simple hashing function which produce an unique numeric 
identifier for each oligo.It is possible to use integer but I use
double to accept a large N.
(See R help:strtoi or just type ?strtoi in console)

*/

inline double  natoi( const std::string A,const int n ){
//Used as input of numeric_split(see next function)
  int np=n-1;
  IntegerVector temp(n,4);
  char toc;
  for( int i=0;i<n;i++){
    toc=A[i];
if(toc=='a') temp[i]= 0;
if(toc=='t') temp[i]= 3;
if(toc=='c') temp[i]= 1;
if(toc=='g') temp[i]= 2;
}
double  result=0;
 for( int i=0;i<n;i++) result+= pow(5,(np-i))*temp[i];
 return result;
}

// [[Rcpp::export]]
inline NumericVector  numeric_split(const std::string A,const int n)
    {
  const int np=n-1;
  const int Lim=A.length()-np;//-- number of subsequence 
  NumericVector res(Lim);     //-- result vector
  double p=0;                 //-- temp result accumulator
  res[0]=p=natoi(A,n);        //-- first result 
  const double  N=pow(5,np);  //-- cont used in recursion formula
  unsigned short int x,y;     //-- temp to allocate forward and backward values
  char fo,ba;                 //--forward and backward chars
  for( int i=1;i<Lim;i++)
     {
       fo=A[i-1];ba=A[i+np];
       x=y=4;                 //-- empty case:n,p,m... 
       if(fo=='a')x=0;        //-- char detector 
       if(fo=='t')x=3;        // 
       if(fo=='g')x=2;        // 
       if(fo=='c')x=1;
       if(ba=='a')y=0;
       if(ba=='t')y=3;
       if(ba=='g')y=2;
       if(ba=='c')y=1;
     res[i]=p= 5*p-5*N*x+y;   //-- scalar formula 
     }
  return res;
  }
/*
R verification :
I use a file named k12.fna in my computer (full E.coli k12 genome ). 
You have just to change k12.fna into your file name.
*/


/***R
x<-read_fasta("k12.fna")  ##change me to your file name 
system.time(charsplit<-SPLIT(x,11))
#user  system elapsed 
# 5.38    0.02    5.39 
system.time(numsplit<-numeric_split(x,11))
#user  system elapsed 
#0.08    0.01    0.09    ## more then 50 time faster
  
  */

Source file example for some sorting algorithms written in C++ via RCPP to R implimentation

#include <Rcpp.h>
# include<algorithm>
using namespace Rcpp;


/*
 * This file contains some algorithms of sorting written in C++ & implimented in RCPP .
 * All these algorithms are not performed as well as sort inner algorithm 
 * used in R except of the sort_stl function .
 *   

 * get pdf: source      
 *                 MACHERKI M.E 
 *                  12/02/2016 
 */


// [[Rcpp::export]]
NumericVector sort_insert(NumericVector A) {
  int i,j;
  double x;
  for(i=1; i<A.size();i++){
    x=A[i];j=i-1;
    while((j>=0) & (A[j]>x)){ //without using swapping
      A[j+1]=A[j];
      j=j-1;
    }
  A[j+1]=x;
  
  }
  return A;
}
 
 // [[Rcpp::export]]
 NumericVector sort_bubble(NumericVector A) { 
   int n=A.size();
   int swap=1;
   double tmp;
   for(;;){
     if(swap==0)break;
     swap=0;
     for(int i=1;i<n;i++){
       if(A[i-1]>A[i]){
         tmp=A[i];A[i]=A[i-1];A[i-1]=tmp;
         swap=1;
       }
     }
   }
   return A;}

// [[Rcpp::export]]
NumericVector sort_selection(NumericVector A) { 
  int i,j;
  double tmp;
  for(j=0;j<A.size();j++){
    int imin=j;
    for(i=j+1;i<A.size();i++){
      if(A[i]<A[imin])imin=i;
    }
  if(imin!=j){
   tmp=A[j];A[j]= A[imin];A[imin]=tmp;
  }
  
  }
  return A;}
  
 

   // [[Rcpp::export]]
 NumericVector gaps( int sz ) { 
   NumericVector result;
   result.push_back(701);
   result.push_back(301);
   result.push_back(132);
   result.push_back(57);
   result.push_back(23);
   result.push_back(10);
   result.push_back(4);
   result.push_back(1);
   return result;}

// [[Rcpp::export]]
NumericVector sort_shell(NumericVector A) { 
  double temp=0;
  int N =A.size();
  int i,j,gap;
  NumericVector tacuda=gaps(0);
  int n=tacuda.size();
  for(int k=n-1;k>0;k--){
    gap=tacuda[k];
  for( i=gap;i<N;i++){
    temp=A[i];
    for( j=i;(j>=gap)&(A[j-gap]>temp);j-=gap){
      A[j]=A[j-gap];
    }
    A[j]=temp;
    
    }
    
  }
  return A;}


  
  
// This function is faster then R sort function !
// [[Rcpp::export]]
NumericVector sort_stl(NumericVector A) { 
  
  NumericVector tmp=clone(A);
  std::sort(tmp.begin(),tmp.end());
  return tmp;
}
  
  
  // [[Rcpp::export]]
 void Quick(NumericVector arr,int left,int right) { 
  int i=left;int j=right;
  double tmp;
  double pivot=arr[(left+right)/2];
  /*partition*/
  while(i<=j){
    while(arr[i]<pivot)
      i++;
    while(arr[j]>pivot)
      j--;
      if(i<=j){
        tmp=arr[i];
        arr[i]=arr[j];
        arr[j]=tmp;
        i++;
        j--;
        
      }
    }
  /*recursion*/
  if(left<j)
    Quick(arr,left,j);
  if(i<right)
    Quick(arr,i,right);
}
    
    
    
    /*
     * This function is really comparable with the inner sort used in R
     * even the code is sample , the algorithm is good.
     */
    
    
    // [[Rcpp::export]]
    NumericVector sort_quick(NumericVector A) {
      NumericVector tmp=clone(A);
    Quick(tmp,0,tmp.size()-1); //just starting from the middle of the vector
      return tmp;
    }
      
  
  
  
  


// R call

/*** R

###not rapid methods, adapted to small sample size
x<-rnorm(100)
sort_insert(x)
sort (x)
sort_selection(x)
sort_bubble(x)
identical(sort_stl(x),sort (x))## the purpose of C++ integration  
identical(sort_quick(x),sort (x))
identical(sort_insert(x),sort (x))
identical(sort_bubble(x),sort (x))
identical(sort_selection(x),sort (x))
identical(sort_shell(x),sort (x))
x<-rnorm(1000000)    ##large numeric vector
system.time(sort(x))
system.time(sort_quick(x)) ## as R sort function
system.time(sort_stl(x)) ## two time faster the R sort!!!    


  */