export to linking to functions and support c++17 parallelization

inst/include/Rfast2/vector.hpp

@@ -11,6 +11,8 @@
 #include "templates.h"
 #include "assertions.hpp"
 
+#include "parallel.h"
+
 namespace Rfast
 {
 
@@ -19,43 +21,46 @@ namespace Rfast
 	using std::remove_if;
 	using std::string;
 
-	template <class Ret, class T>
-	Ret Quantile(T x, colvec &probs)
+	template <class Ret, class T, class F>
+	Ret Quantile(T x, F &probs, const bool parallel = false)
 	{
 		Assertion::is_iterable<T>::check_concept();
 		Assertion::has_size<T>::check_concept();
-
-		Ret f(probs.n_elem);
+		Assertion::is_iterable<F>::check_concept();
+		Assertion::has_size<F>::check_concept();
+
+		const unsigned int nprobs = probs.size();
+
+		Ret f(nprobs);
 
-		if (probs.size() > std::log2(x.size()))
+		if (nprobs > std::log2(x.size()))
 		{ // k > log2(n) tote allazo algorithmo
 			const int mxelem = (x.n_elem - 1) * (*max_element(probs.begin(), probs.end())) + 1;
 			std::nth_element(x.begin(), x.begin() + mxelem, x.end());
-			std::sort(x.begin(), x.end());
-			for (unsigned int i = 0; i < probs.n_elem; ++i)
+			Rfast::sort(x.begin(), x.end(), parallel);
+			for (unsigned int i = 0; i < nprobs; ++i)
 			{
 				double h = (x.n_elem - 1) * probs[i] + 1;
 				int hf = h;
 				auto a = x[hf - 1];
-				auto b = x[hf];
-				f[i] = a + (h - hf) * (b - a);
+				f[i] = a + (h - hf) * (x[hf] - a);
 			}
 		}
 		else
 		{
-			for (unsigned int i = 0; i < probs.n_elem; ++i)
+			for (unsigned int i = 0; i < nprobs; ++i)
 			{
-				double h = (x.n_elem - 1) * probs[i] + 1;
+				double h = (x.size() - 1) * probs[i] + 1;
 				int hf = h;
 				double a, b;
 				if (probs[i] > 0.5)
 				{
-					a = nth_simple<T>(x, hf - 1, false);
+					a = nth_simple<T>(x, hf - 1, false, parallel);
 					b = *min_element(x.begin() + hf, x.end());
 				}
 				else
 				{
-					b = nth_simple<T>(x, hf, false);
+					b = nth_simple<T>(x, hf, false, parallel);
 					a = *max_element(x.begin(), x.begin() + hf);
 				}
 				f[i] = a + (h - hf) * (b - a);
@@ -65,27 +70,19 @@ namespace Rfast
 		return f;
 	}
 
-
-
 	template <class T>
-	double TrimMean(T x, const double a)
+	double TrimMean(T x, const double a = 0.5, const bool parallel = false)
 	{
 		Assertion::is_iterable<T>::check_concept();
 		Assertion::has_size<T>::check_concept();
 
 		const int n = x.size();
 		int b1 = a * n;
 		int b11 = std::ceil(b1);
-		if (b1 - b11 == 0)
-		{
-			b1 = b11 + 1;
-		}
-		else
-		{
-			b1 = b11;
-		}
-		const double a1 = nth_simple<T>(x, b1, false);
-		const double a2 = nth_simple<T>(x, n - b1 + 1, false);
+		b1 = (b1 == b11) ? b11 + 1 : b11;
+
+		const double a1 = nth_simple<T>(x, b1, false, parallel);
+		const double a2 = nth_simple<T>(x, n - b1 + 1, false, parallel);
 		double s = 0;
 		int p = 0;
 		for (auto xx : x)

src/utilities.cpp

@@ -144,28 +144,30 @@ END_RCPP
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-RcppExport SEXP Rfast2_Quantile(SEXP xSEXP,SEXP ProbsSEXP){
+RcppExport SEXP Rfast2_Quantile(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP){
 BEGIN_RCPP
     RObject __result;
     RNGScope __rngScope;
+    traits::input_parameter< const bool >::type parallel(parallelSEXP);
     NumericVector x(xSEXP);
     NumericVector Probs(ProbsSEXP);
 
     colvec probs(Probs.begin(),Probs.size(),false);
-    __result = Rfast::Quantile<NumericVector,colvec>(x,probs);
+    __result = Rfast::Quantile<NumericVector,colvec>(x,probs,parallel);
     return __result;
 END_RCPP
 }
 
 /***********************************************************************************/
 
-RcppExport SEXP Rfast2_trimmean(SEXP xSEXP,SEXP aSEXP){
+RcppExport SEXP Rfast2_trimmean(SEXP xSEXP,SEXP aSEXP,SEXP parallelSEXP){
 BEGIN_RCPP
     RObject __result;
     RNGScope __rngScope;
     traits::input_parameter< colvec  >::type x(xSEXP);
     traits::input_parameter< const double >::type a(aSEXP);
-    __result = Rfast::TrimMean<colvec>(x,a);
+    traits::input_parameter< const bool >::type parallel(parallelSEXP);
+    __result = Rfast::TrimMean<colvec>(x,a,parallel);
     return __result;
 END_RCPP
 }