diff --git a/.Rbuildignore b/.Rbuildignore index f7f0ab7..3af075d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -3,4 +3,6 @@ ^docs ^.todo ^.github -^.html \ No newline at end of file +^.html +^.gitignore +^.vscode \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..7aeb440 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,81 @@ +{ + "files.associations": { + "random": "cpp", + "array": "cpp", + "chrono": "cpp", + "deque": "cpp", + "format": "cpp", + "forward_list": "cpp", + "initializer_list": "cpp", + "list": "cpp", + "vector": "cpp", + "xhash": "cpp", + "xstring": "cpp", + "xtree": "cpp", + "xutility": "cpp", + "algorithm": "cpp", + "utility": "cpp", + "atomic": "cpp", + "bit": "cpp", + "cctype": "cpp", + "charconv": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "compare": "cpp", + "concepts": "cpp", + "cstddef": "cpp", + "cstdint": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "exception": "cpp", + "fstream": "cpp", + "iomanip": "cpp", + "ios": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "iterator": "cpp", + "limits": "cpp", + "locale": "cpp", + "map": "cpp", + "memory": "cpp", + "mutex": "cpp", + "new": "cpp", + "optional": "cpp", + "ostream": "cpp", + "ratio": "cpp", + "sstream": "cpp", + "stack": "cpp", + "stdexcept": "cpp", + "stop_token": "cpp", + "streambuf": "cpp", + "string": "cpp", + "system_error": "cpp", + "thread": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "typeinfo": "cpp", + "unordered_map": "cpp", + "xfacet": "cpp", + "xiosbase": "cpp", + "xlocale": "cpp", + "xlocbuf": "cpp", + "xlocinfo": "cpp", + "xlocmes": "cpp", + "xlocmon": "cpp", + "xlocnum": "cpp", + "xloctime": "cpp", + "xmemory": "cpp", + "xstddef": "cpp", + "xtr1common": "cpp", + "functional": "cpp", + "queue": "cpp", + "*.rmd": "markdown", + "unordered_set": "cpp", + "execution": "cpp", + "numeric": "cpp" + } +} \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 0c403e9..94ffd5b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,13 @@ Package: Rfast2 Type: Package Title: A Collection of Efficient and Extremely Fast R Functions II -Version: 0.1.4 -Date: 2023-02-14 +Version: 0.1.5.1 +Date: 2023-07-21 Author: Manos Papadakis, Michail Tsagris, Stefanos Fafalios and Marios Dimitriadis. Maintainer: Manos Papadakis -Depends: R (>= 3.5.0), Rcpp (>= 0.12.3) -LinkingTo: Rcpp (>= 0.12.3), RcppArmadillo +Depends: R (>= 3.5.0), Rcpp (>= 0.12.3), RcppParallel +LinkingTo: Rcpp (>= 0.12.3), RcppArmadillo, RcppParallel +SystemRequirements: C++17 Imports: Rfast, RANN BugReports: https://github.com/RfastOfficial/Rfast2/issues URL: https://github.com/RfastOfficial/Rfast2 diff --git a/NAMESPACE b/NAMESPACE index 275ee77..c608b7d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,7 @@ importFrom("Rfast", vm.mle, colvm.mle, colmeans, colsums, eachrow,group, colVars importFrom("RANN", nn2) importFrom(Rcpp, evalCpp) +importFrom(RcppParallel, RcppParallelLibs) export(.onAttach,Intersect,Merge,Quantile,Runif,Sample,Sample.int,add.term,benchmark,bernoulli.nb,bernoullinb.pred,beta.nb,betanb.pred,bic.regs,big.knn,bigknn.cv,binom.reg,boot.hotel2,boot.james,boot.student2,boot.ttest1,cauchy.nb,cauchy0.mle,cauchy0.mle,cauchynb.pred,cbern.mle,censpois.mle,censpois.mle,censweib.reg,censweibull.mle,censweibull.mle,circ.cor1,circ.cors1,cls,cluster.lm,col.waldpoisrat,colGroup,colQuantile,colTrimMean,colaccs,colbeta.mle,colborel.mle,colcauchy.mle,colfbscores,colfmis,colfscores,colhalfnorm.mle,coljack.means,collogitnorm.mle,collognorm.mle,colmaes,colmeansvars,colmses,colordinal.mle,colpinar1,colpkl,colprecs,colsens,colspecs,colspml.mle,colukl,colunitweibull.mle,colwlsmeta,cor_test,covar,covdist,covequal,covlikel,covmtest,covrob.lm,dcora,den.contours,depth.mahala,diffic,discrim,embed.circaov,empirical.entropy,fbed.reg,fe.lmfit,fedhc.skel,fipois.reg,fisher.da,fp,gammapois.mle,gammapois.mle,gammareg,gammaregs,gee.reg,gnormal0.mle,gumbel.reg,halfcauchy.mle,halfcauchy.mle,hcf.circaov,hellinger.countreg,het.circaov,het.lmfit,hp.reg,is.lower.tri,is.skew.symmetric,is.upper.tri,jack.mean,km,kumar.mle,kumar.mle,laplace.nb,laplacenb.pred,leverage,lm.bsreg,lm.drop1,lm.parboot,logiquant.regs,logitnorm.nb,logitnormnb.pred,lr.circaov,lud,mci,mle.lda,mmhc.skel,mmpc,mmpc2,moranI,multinom.reg,multinomreg.cv,multispml.mle,multivm.mle,mv.score.betaregs,mv.score.expregs,mv.score.gammaregs,mv.score.glms,mv.score.invgaussregs,mv.score.weibregs,nb.cv,negbin.reg,negbin.regs,normlog.nb,normlognb.pred,omp2,overdispreg.test,pc.sel,pca,pcr,perm.ttest2,pinar1,pooled.colVars,powerlaw.mle,powerlaw.mle,prophelling.reg,propols.reg,purka.mle,purka.mle,rbeta1,refmeta,reg.mle.lda,regmlelda.cv,riag,rm.hotel,rowQuantile,rowTrimMean,rowjack.means,sclr,score.zipregs,simplex.mle,simplex.mle,sp.logiregs,sp.mle,spml.nb,spmlnb.pred,stud.ttests,tobit.reg,trim.mean,trunccauchy.mle,trunccauchy.mle,truncexpmle,truncexpmle,unitweibull.mle,vm.nb,vmnb.pred,wald.poisrat,walter.ci,weib.regs,weibull.nb,weibullnb.pred,welch.tests,wild.boot,wlsmeta,zigamma.mle,zigamma.mle,zigamma.reg,zil.mle,zil.mle,ziweibull.mle,ziweibull.mle,ztp.reg) diff --git a/NEWS.md b/NEWS.md index 204ce02..ab0856e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,219 +1,247 @@ -

_Rfast2_

+

-### **Version 0.1.4** -#### Date release: **01/11/2022** -*** -> **New** +*Rfast2* + +

+ +### **Version 0.1.5 - Theseus** + +------------------------------------------------------------------------ + +>

**Improved**

(***by speed, correctness or options***) > -> | Function | What's new! | -> | :----------------: | :------------------------------------------------: | -> | covrob.lm | Linear regression with robust covariance matrix | +> | Function | What's new! | +> |:--------:|:-----------------------------------------------------------------------------------:| +> | TrimMean | Add option for parallelism. Supported only where C++ execution policy is supported. | +> | Quantile | Add option for parallelism. Supported only where C++ execution policy is supported. | > +>

**LinkingTo**

(***by speed, correctness or options***) > -> **Improved** (_**by speed, correctness or options**_) +> | Function/Structure | What's new! | +> |:------------------:|:--------------------:| +> | Quantile | Exported for linking | +> | rowQuantile | Exported for linking | +> | colQuantile | Exported for linking | +> | TrimMean | Exported for linking | +> | rowTrimMean | Exported for linking | +> | colTrimMean | Exported for linking | +> | colTrimMean | Exported for linking | > -> | Function | What's new! | -> | :-----------------: | :---------------------------------------: | -> | | | +> ### **Comments** +> +> ------------------------------------------------------------------------ +> +> - All the exported functions for linking to mechanism are in namespace `Rfast`. +> - We have added the header file `parallel.h` which automatically includes the new parallelism policy rules of C++17. If not supported by your system, non-parallel algorithms are automatically used. + +

+ +### **Version 0.1.4 - Heracles** -
-
+------------------------------------------------------------------------ + +>

**New**

+> +> | Function | What's new! | +> |:---------:|:-----------------------------------------------:| +> | covrob.lm | Linear regression with robust covariance matrix | + +

### **Version 0.1.3** -#### Date release: **22/03/2022** -*** -> **New** -> -> | Function | What's new! | -> | :----------------: | :------------------------------------------------: | -> | Runif | Like R's runif but faster. | -> | Sample | Like R's Sample but faster. | -> | Sample.int | Like R's Sample.int but faster. | -> | colaccs | Column-wise accuracies | -> | colsens | Column-wise sensitivities | -> | colspecs | Column-wise specificities | -> | colprecs | Column-wise precisions | -> | colfscores | Column-wise F-scores | -> | colfbscores | Column-wise F-beta-scores | -> | colfmis | Column-wise Fowlkes–Mallows index | -> | colfmses | Column-wise MSEs | -> | colmaes | Column-wise MAEs | -> | colpkl | Column-wise Kullback-Leibler divergence for percentages | -> | colukl | Column-wise Kullback-Leibler divergence for non-negative or non-negative values | -> | pinar1 | Poisson INAR(1) model estimation | -> | colpinar1 | Column-wise Poisson INAR(1) model estimation | -> -> -> **Improved** (_**by speed, correctness or options**_) -> -> | Function | What's new! | -> | :-----------------: | :-----------------------------------------: | -> | Quantile, rowQuantile | Optimize algorithm | -> | colQuantile | Optimize algorithm and new method for data.frames | -> | colTrimMean | Optimize algorithm and new method for data.frames | - -
-
+ +------------------------------------------------------------------------ + +>

**New**

+> +> | Function | What's new! | +> |:-----------:|:-------------------------------------------------------------------------------:| +> | Runif | Like R's runif but faster. | +> | Sample | Like R's Sample but faster. | +> | Sample.int | Like R's Sample.int but faster. | +> | colaccs | Column-wise accuracies | +> | colsens | Column-wise sensitivities | +> | colspecs | Column-wise specificities | +> | colprecs | Column-wise precisions | +> | colfscores | Column-wise F-scores | +> | colfbscores | Column-wise F-beta-scores | +> | colfmis | Column-wise Fowlkes--Mallows index | +> | colfmses | Column-wise MSEs | +> | colmaes | Column-wise MAEs | +> | colpkl | Column-wise Kullback-Leibler divergence for percentages | +> | colukl | Column-wise Kullback-Leibler divergence for non-negative or non-negative values | +> | pinar1 | Poisson INAR(1) model estimation | +> | colpinar1 | Column-wise Poisson INAR(1) model estimation | +> +>

**Improved**

(***by speed, correctness or options***) +> +> | Function | What's new! | +> |:---------------------:|:-------------------------------------------------:| +> | Quantile, rowQuantile | Optimize algorithm | +> | colQuantile | Optimize algorithm and new method for data.frames | +> | colTrimMean | Optimize algorithm and new method for data.frames | + +

### **Version 0.1.1** -#### Date release: **10/10/2021** -*** -> **New** + +------------------------------------------------------------------------ + +>

**New**

> -> | Function | What's new! | -> | :----------------: | :------------------------------------------------: | -> | mmhc.skel | Skeleton of MMHC Bayesian network learning algorithm | -> | fedhc.skel | Skeleton of FEDHC Bayesian network learning algorithm | -> | fe.lmfit | Fixed effects linear regression for panel data | +> | Function | What's new! | +> |:----------:|:-----------------------------------------------------:| +> | mmhc.skel | Skeleton of MMHC Bayesian network learning algorithm | +> | fedhc.skel | Skeleton of FEDHC Bayesian network learning algorithm | +> | fe.lmfit | Fixed effects linear regression for panel data | -
-
+

### **Version 0.0.8** -#### Date release: **16/12/2019** -*** -> **New** -> -> | Function | What's new! | -> | :----------------: | :------------------------------------------------: | -> | gammareg | Gamma regression | -> | gammaregs | Many gamma regressions | -> | cor_test | Hypothesis testing for the correlation coefficient | -> | cor_test | Hypothesis testing for the correlation coefficient | -> | hcf.circaov | ANOVA for circular data | -> | het.circaov | ANOVA for circular data | -> | lr.circaov | ANOVA for circular data | -> | multivm.mle | Fitting many von Mises distributions | -> | multispml.mle | Fitting many von SPML distributions | -> | depth.mahala | Mahalanobis depth | -> | perm.ttest2 | Permutation based 2 sample t-test | -> | lm.parboot | Parametric bootstrap for linear models | -> | weibull.nb | Weibull Naive Bayes (NB) classifier | -> | weibullnb.pred | Prediction using Weibull NB classifier | -> | normlog.nb | Gaussian(log-link) NB classifier | -> | normlognb.pred | Prediction using Gaussian(log-link) NB classifier | -> | laplace.nb | Laplace NB classifier | -> | laplacenb.pred | Prediction using Laplace NB classifier | -> | vm.nb | von Mises NB classifier | -> | vmnb.pred | Prediction using von Mises NB classifier | -> | spml.nb | SPML NB classifier | -> | spml.pred | Prediction using SPML NB classifier | - -
-
+ +------------------------------------------------------------------------ + +>

**New**

+> +> | Function | What's new! | +> |:--------------:|:--------------------------------------------------:| +> | gammareg | Gamma regression | +> | gammaregs | Many gamma regressions | +> | cor_test | Hypothesis testing for the correlation coefficient | +> | cor_test | Hypothesis testing for the correlation coefficient | +> | hcf.circaov | ANOVA for circular data | +> | het.circaov | ANOVA for circular data | +> | lr.circaov | ANOVA for circular data | +> | multivm.mle | Fitting many von Mises distributions | +> | multispml.mle | Fitting many von SPML distributions | +> | depth.mahala | Mahalanobis depth | +> | perm.ttest2 | Permutation based 2 sample t-test | +> | lm.parboot | Parametric bootstrap for linear models | +> | weibull.nb | Weibull Naive Bayes (NB) classifier | +> | weibullnb.pred | Prediction using Weibull NB classifier | +> | normlog.nb | Gaussian(log-link) NB classifier | +> | normlognb.pred | Prediction using Gaussian(log-link) NB classifier | +> | laplace.nb | Laplace NB classifier | +> | laplacenb.pred | Prediction using Laplace NB classifier | +> | vm.nb | von Mises NB classifier | +> | vmnb.pred | Prediction using von Mises NB classifier | +> | spml.nb | SPML NB classifier | +> | spml.pred | Prediction using SPML NB classifier | + +

### **Version 0.0.4** -#### Date release: **28/07/2019** -*** -> **Improved** (_**by speed, correctness or options**_) -> -> | Function | What's new! | -> | :-----------------: | :-----------------------------------------: | -> | bic.regs | Addition SPML and Weibull regression | -> | pc.sel | Time improvement | -> | welch.tests | Time improevement | -> -> -> **New** -> -> | Function | What's new! | -> | :-----------------: | :---------------------------------------------------: | -> | cls | Constrained least squares regression | -> | cluster.lm | Linear regression with clustered data | -> | colborel.mle | Column-wise MLE of the Borel distribution | -> | collogitnorm.mle | Column-wise MLE of the logistic normal distribution | -> | collognorm.mle | Column-wise MLE of the log-normal distribution | -> | cospml.mle | Column-wise MLE of the SPML distribution | -> | empirical.entropy | Empirical entropy estimation with continuous data | -> | fbed.reg | FBED variable selection algorithm | -> | gumbel.reg | Gumbel regression | -> | moranI | Moran's I measure of spatial autocorrelation | -> | multinom.reg | Multinomial regression | -> | negbin.reg | Negative binomial regression | -> | pca | Principal component analysis | -> | refmeta | Random effects meta-analysis estimation | -> | simplex.mle | MLE of the simplex distribution | -> | weib.regs | Many simple weibull regressions | -> | ztp.reg | zero truncated Poisson regression | -> | mmpc2 | MMPC variable selection algorithm | -> | is.skew.symmetric | Checking if the given matrix is skew-symmetric. | - -
-
+ +------------------------------------------------------------------------ + +>

**Improved**

(***by speed, correctness or options***) +> +> | Function | What's new! | +> |:-----------:|:------------------------------------:| +> | bic.regs | Addition SPML and Weibull regression | +> | pc.sel | Time improvement | +> | welch.tests | Time improevement | +> +>

**New**

+> +> | Function | What's new! | +> |:-----------------:|:---------------------------------------------------:| +> | cls | Constrained least squares regression | +> | cluster.lm | Linear regression with clustered data | +> | colborel.mle | Column-wise MLE of the Borel distribution | +> | collogitnorm.mle | Column-wise MLE of the logistic normal distribution | +> | collognorm.mle | Column-wise MLE of the log-normal distribution | +> | cospml.mle | Column-wise MLE of the SPML distribution | +> | empirical.entropy | Empirical entropy estimation with continuous data | +> | fbed.reg | FBED variable selection algorithm | +> | gumbel.reg | Gumbel regression | +> | moranI | Moran's I measure of spatial autocorrelation | +> | multinom.reg | Multinomial regression | +> | negbin.reg | Negative binomial regression | +> | pca | Principal component analysis | +> | refmeta | Random effects meta-analysis estimation | +> | simplex.mle | MLE of the simplex distribution | +> | weib.regs | Many simple weibull regressions | +> | ztp.reg | zero truncated Poisson regression | +> | mmpc2 | MMPC variable selection algorithm | +> | is.skew.symmetric | Checking if the given matrix is skew-symmetric. | + +

### **Version 0.0.3** -#### Date release: **08/05/2019** -*** -> **Improved** (_**by speed, correctness or options**_) + +------------------------------------------------------------------------ + +>

**Improved**

(***by speed, correctness or options***) > -> | Function | What's new! | -> | :-----------------: | :-----------------------------------------: | -> | benchmark | fix a bug in printing names. | +> | Function | What's new! | +> |:---------:|:----------------------------:| +> | benchmark | fix a bug in printing names. | -
-
+

### **Version 0.0.2** -#### Date release: **08/05/2019** -*** -> **New** -> -> | Function | What's new! | -> | :-----------------: | :----------------------------------------------------: | -> | gee.reg | Gereneralised estimating equations Gaussian regression | -> | halfcauchy.mle | MLE of the half Cauchy distribution | -> | kumar.mle | MLE of the Kumaraswamy distribution | -> | powerlaw.mle | MLE of the power law distribution | -> | reg.mle.lda | Regularised maximum likleihood linear discriminant analysis | -> | walter.ci | Confidence interval for the relative risk using Walter's method | -> | welch.tests | Many Welch t-tests | - -
-
+ +------------------------------------------------------------------------ + +>

**New**

+> +> | Function | What's new! | +> |:--------------:|:---------------------------------------------------------------:| +> | gee.reg | Gereneralised estimating equations Gaussian regression | +> | halfcauchy.mle | MLE of the half Cauchy distribution | +> | kumar.mle | MLE of the Kumaraswamy distribution | +> | powerlaw.mle | MLE of the power law distribution | +> | reg.mle.lda | Regularised maximum likleihood linear discriminant analysis | +> | walter.ci | Confidence interval for the relative risk using Walter's method | +> | welch.tests | Many Welch t-tests | + +

### **Version 0.0.1** -#### Date release: **19/03/2019** -*** -> **New** -> -> | Function | What's new! | -> | :-----------------: | :-------------------------------------------------------------: | -> | add.term | Add many single terms to a model. | -> | bic.regs | BIC of many univariate regressions. | -> | censpois.mle | MLE of the left censored Poisson distribution. | -> | cesweibull.mle | MLE of the censored Weibull distribution. | -> | circ.cor1 | Circurlar correlations between two circular variables. | -> | circ.cors1 | Circurlar correlations between a circular and many circular variables. | -> | colmeansvars | Column-wise means and variances of a matrix. | -> | covar | Covariance betweeen a vector and a matrix. | -> | diffic | Difficulty of items (psychometric theory). | -> | discrim | Discrimination of items (psychometric theory). | -> | gammapois.mle | MLE of the gamma-Poisson distribution. | -> | km | Kaplan-Meier estimate of a survival function. | -> | logiquant.regs | Many simple quantile regressions using logistic regressions. | -> | mle.lda | Maximum likelihood linear discriminant analysis. | -> | pc.sel | Variable selection using the PC-simple algorithm. | -> | pooled.colVars | Column-wise pooled variances across groups. | -> | purka.mle | MLE of the Purkayastha distribution. | -> | sp.logiregs | Many approximate simple logistic regressionss. | -> | trunccauchy.mle | MLE of the truncated Cauchy distribution. | -> | truncexpmle | MLE of the truncated exponential distribution. | -> | wald.poisrat | Wald confidence interval for the ratio of two Poisson variables. | -> | col.waldpoisrat | Column-wise Wald confidence interval for the ratio of two Poisson variables. | -> | welch.tests | Many Welch tests. | -> | zigamma.mle | MLE of the zero inflated Gamma distribution. | -> | ziweibull.mle | MLE of the zero inflated Weibull distribution. | -> | zil.mle | MLE of the zero inflated logistic normal distribution. | -> | Intersect | Intersect as R's. | -> | is.lower.tri | Check if a matrix is lower triangular. | -> | is.upper.tri | Check if a matrix is upper triangular. | -> | lud | Split a matrix to a lower,upper matrix and diagonal vector.| -> | Merge | Merge 2 sorted vectors to a sorted vector. | -> | benchmark | Measure code's execution time. | -> | colGroup | Apply Rfast's gorup function to each column with some restrictions. | -> | Quantile | Quantile(s) of a vector. | -> | colQuantile | Column-wise quantile(s) of a matrix. | -> | rowQuantile | Row-wise quantile(s) of a matrix. | -> | trim.mean | Trimmed mean of a vector. | -> | colTrimMean | Column-wise trimmed mean of a matrix. | -> | rowTrimMean | Row-wise trimmed mean of a matrix. | + +------------------------------------------------------------------------ + +>

**New**

+> +> | Function | What's new! | +> |:---------------:|:----------------------------------------------------------------------------:| +> | add.term | Add many single terms to a model. | +> | bic.regs | BIC of many univariate regressions. | +> | censpois.mle | MLE of the left censored Poisson distribution. | +> | cesweibull.mle | MLE of the censored Weibull distribution. | +> | circ.cor1 | Circurlar correlations between two circular variables. | +> | circ.cors1 | Circurlar correlations between a circular and many circular variables. | +> | colmeansvars | Column-wise means and variances of a matrix. | +> | covar | Covariance betweeen a vector and a matrix. | +> | diffic | Difficulty of items (psychometric theory). | +> | discrim | Discrimination of items (psychometric theory). | +> | gammapois.mle | MLE of the gamma-Poisson distribution. | +> | km | Kaplan-Meier estimate of a survival function. | +> | logiquant.regs | Many simple quantile regressions using logistic regressions. | +> | mle.lda | Maximum likelihood linear discriminant analysis. | +> | pc.sel | Variable selection using the PC-simple algorithm. | +> | pooled.colVars | Column-wise pooled variances across groups. | +> | purka.mle | MLE of the Purkayastha distribution. | +> | sp.logiregs | Many approximate simple logistic regressionss. | +> | trunccauchy.mle | MLE of the truncated Cauchy distribution. | +> | truncexpmle | MLE of the truncated exponential distribution. | +> | wald.poisrat | Wald confidence interval for the ratio of two Poisson variables. | +> | col.waldpoisrat | Column-wise Wald confidence interval for the ratio of two Poisson variables. | +> | welch.tests | Many Welch tests. | +> | zigamma.mle | MLE of the zero inflated Gamma distribution. | +> | ziweibull.mle | MLE of the zero inflated Weibull distribution. | +> | zil.mle | MLE of the zero inflated logistic normal distribution. | +> | Intersect | Intersect as R's. | +> | is.lower.tri | Check if a matrix is lower triangular. | +> | is.upper.tri | Check if a matrix is upper triangular. | +> | lud | Split a matrix to a lower,upper matrix and diagonal vector. | +> | Merge | Merge 2 sorted vectors to a sorted vector. | +> | benchmark | Measure code's execution time. | +> | colGroup | Apply Rfast's gorup function to each column with some restrictions. | +> | Quantile | Quantile(s) of a vector. | +> | colQuantile | Column-wise quantile(s) of a matrix. | +> | rowQuantile | Row-wise quantile(s) of a matrix. | +> | trim.mean | Trimmed mean of a vector. | +> | colTrimMean | Column-wise trimmed mean of a matrix. | +> | rowTrimMean | Row-wise trimmed mean of a matrix. | diff --git a/R/Quantile.R b/R/Quantile.R index 923b6c4..abba181 100644 --- a/R/Quantile.R +++ b/R/Quantile.R @@ -1,22 +1,22 @@ -#[export] -Quantile <- function(x, probs) { - .Call(Rfast2_Quantile, x, probs) +# [export] +Quantile <- function(x, probs, parallel = FALSE) { + .Call(Rfast2_Quantile, x, probs, parallel) } -#[export] -rowQuantile <- function(x, probs, parallel = FALSE) { - .Call(Rfast2_rowQuantile, x, probs, parallel) +# [export] +rowQuantile <- function(x, probs, parallel = FALSE, cores = 0) { + .Call(Rfast2_rowQuantile, x, probs, parallel, cores) } -#[export s3] -colQuantile.data.frame<-function(x, probs, parallel = FALSE) { - .Call(Rfast2_colQuantile, x, probs, parallel) +# [export s3] +colQuantile.data.frame <- function(x, probs, parallel = FALSE, cores = 0) { + .Call(Rfast2_colQuantile, x, probs, parallel, cores) } -#[export] -colQuantile<-function(x, probs, parallel = FALSE) { - UseMethod("colQuantile") +# [export] +colQuantile <- function(x, probs, parallel = FALSE, cores = 0) { + UseMethod("colQuantile") } -#[export s3] -colQuantile.matrix<-function(x, probs, parallel = FALSE) { - .Call(Rfast2_colQuantile, x, probs, parallel) +# [export s3] +colQuantile.matrix <- function(x, probs, parallel = FALSE, cores = 0) { + .Call(Rfast2_colQuantile, x, probs, parallel, cores) } diff --git a/R/onAttach.R b/R/onAttach.R index 7b5edc7..157d0cc 100644 --- a/R/onAttach.R +++ b/R/onAttach.R @@ -1,5 +1,5 @@ -#[export] +#[export special] .onAttach <- function(lib, pkg) { version <- read.dcf(file.path(lib, pkg, "DESCRIPTION"),"Version") packageStartupMessage(paste("\nRfast2: ",version)) diff --git a/R/trim.mean.R b/R/trim.mean.R index 4ed4e8e..74ae54a 100644 --- a/R/trim.mean.R +++ b/R/trim.mean.R @@ -1,24 +1,24 @@ #[export] -trim.mean <- function(x, a = 0.05) { - .Call(Rfast2_trimmean, x, a) +trim.mean <- function(x, a = 0.05, parallel = FALSE) { + .Call(Rfast2_trimmean, x, a, parallel) } #[export] -rowTrimMean <- function(x, a = 0.05, parallel = FALSE) { - .Call(Rfast2_rowTrimMean, x, a, parallel) +rowTrimMean <- function(x, a = 0.05, parallel = FALSE, cores = 0) { + .Call(Rfast2_rowTrimMean, x, a, parallel, cores) } #[export] -colTrimMean <- function(x, a = 0.05, parallel = FALSE) { - UseMethod("colTrimMean") +colTrimMean <- function(x, a = 0.05, parallel = FALSE, cores = 0) { + UseMethod("colTrimMean") } #[export s3] -colTrimMean.matrix <- function(x, a = 0.05, parallel = FALSE) { - .Call(Rfast2_colTrimMean, x, a, parallel) +colTrimMean.matrix <- function(x, a = 0.05, parallel = FALSE, cores = 0) { + .Call(Rfast2_colTrimMean, x, a, parallel, cores) } #[export s3] -colTrimMean.data.frame <- function(x, a = 0.05, parallel = FALSE) { - .Call(Rfast2_colTrimMean, x, a, parallel) +colTrimMean.data.frame <- function(x, a = 0.05, parallel = FALSE, cores = 0) { + .Call(Rfast2_colTrimMean, x, a, parallel, cores) } diff --git a/R/univariate.mle.R b/R/univariate.mle.R index a1fb528..6e3d70a 100644 --- a/R/univariate.mle.R +++ b/R/univariate.mle.R @@ -442,52 +442,6 @@ purka.mle <- function(x, tol = 1e-07) { } -#[export] -simplex.mle <- function(x, tol = 1e-07) { - n <- length(x) - xx <- x * (1 - x) - simplexfun <- function(m, xx, x) sum( (x - m)^2 /xx ) / ( m^2 * (1 - m)^2 ) - mod <- optimise(simplexfun, c(0, 1), xx = xx, x = x, tol = tol) - s <- sqrt( mod$objective/n ) - param <- c( mod$minimum, s) - names(param) <- c("mean", "sigma") - list(param = param, loglik = -0.5 * n * log(2 * pi) - 1.5 * sum( log(xx) ) - n * log(s) - n/2 ) -} - -# simplex.mle <- function (x, tol = 1e-09) { -# n <- length(x) -# xx <- x * (1 - x) -# a <- min(x) -# b <- max(x) -# ratio <- 2/(sqrt(5) + 1) -# m1 <- b - ratio * (b - a) -# m2 <- a + ratio * (b - a) -# f1 <- - sum( (x - m1)^2 /xx ) / ( m1^2 * (1 - m1)^2 ) -# f2 <- - sum( (x - m2)^2 /xx ) / ( m2^2 * (1 - m2)^2 ) -# -# while ( abs(f2 - f1) > tol ) { -# if (f2 < f1) { -# b <- m2 -# m2 <- m1 -# f2 <- f1 -# m1 <- b - ratio * (b - a) -# f1 <- - sum( (x - m1)^2 /xx ) / ( m1^2 * (1 - m1)^2 ) -# } else { -# a <- m1 -# m1 <- m2 -# f1 <- f2 -# m2 <- a + ratio * (b - a) -# f2 <- - sum( (x - m2)^2 /xx ) / ( m2^2 * (1 - m2)^2 ) -# } -# } -# -# m <- 0.5 * (m1 + m2) -# s <- sqrt( - f2/n ) -# param <- c( m, s) -# names(param) <- c("mean", "sigma") -# list(param = param, loglik = -0.5 * n * log(2 * pi) - 1.5 * sum( log(xx) ) - n * log(s) - n/2 ) -# } - #[export] trunccauchy.mle <- function (x, a, b, tol = 1e-07) { diff --git a/inst/include/Rfast2.h b/inst/include/Rfast2.h new file mode 100644 index 0000000..aa910ba --- /dev/null +++ b/inst/include/Rfast2.h @@ -0,0 +1,8 @@ + +#ifndef RFAST2_H +#define RFAST2_H + +#include "Rfast2/matrix.hpp" +#include "Rfast2/vector.hpp" + +#endif \ No newline at end of file diff --git a/inst/include/Rfast2/assertions.hpp b/inst/include/Rfast2/assertions.hpp new file mode 100644 index 0000000..7c8f7c0 --- /dev/null +++ b/inst/include/Rfast2/assertions.hpp @@ -0,0 +1,46 @@ + + +#ifndef ASSERTIONS_H +#define ASSERTIONS_H + +namespace Assertion +{ + +#include +#include + + using namespace std; + + template + struct is_iterable : false_type + { + constexpr static void check_concept() + { + static_assert(!is_same::value, "The template class must be iterable."); + } + }; + + template + struct is_iterable().begin()), decltype(declval().end())>> : true_type + { + constexpr static void check_concept() {} + }; + + template + struct has_size : false_type + { + constexpr static void check_concept() + { + static_assert(!is_same::value, "The template class must provide a function named size."); + } + }; + + template + struct has_size().size())>> : true_type + { + constexpr static void check_concept() {} + }; +} + +#endif \ No newline at end of file diff --git a/inst/include/Rfast2/helpers.hpp b/inst/include/Rfast2/helpers.hpp new file mode 100644 index 0000000..0500b8d --- /dev/null +++ b/inst/include/Rfast2/helpers.hpp @@ -0,0 +1,259 @@ + +#ifndef HELPERS_HPP +#define HELPERS_HPP + +#include +#include +#ifdef _OPENMP +#include +#endif + +using namespace std; +using namespace arma; +using namespace Rcpp; + +inline unsigned int get_num_of_threads() +{ +#ifdef _OPENMP + return omp_get_max_threads(); +#else + return 0; +#endif +} + +template Function, class FF> +typename T::value_type parallelSingleIteratorWithoutCopy(FF s) +{ + T y; +#pragma omp critical + { + HELPER yy(s); + y = T(yy.begin(), yy.size(), false); + } + return *Function(y.begin(), y.end()); +} + +template Function, class FF> +typename T::value_type singleIteratorWithoutCopy(FF c) +{ + Helper h(c); + T y(h.begin(), h.size(), false); + return *Function(y.begin(), y.end()); +} + +template Function, class FF> +void singleIteratorWithoutCopy(List &f, const int i, FF c) +{ + T y(c); + f[i] = T::create(*Function(y.begin(), y.end())); +} + +template Function, class FF> +void singleIteratorWithoutCopy(mat &f, const int i, FF c) +{ + Helper h(c); + T y(h.begin(), h.size(), false); + f[i] = *Function(y.begin(), y.end()); +} + +template Function, class FF> +void setResult(List &f, const int i, FF c) +{ + T y = clone(as(c)); + Function(y.begin(), y.end()); + f[i] = T(y.begin(), y.end()); +} + +template ::type, + typename std::remove_reference::type>> + Function, + class FF, + class F> +void setResult(mat &f, const int i, FF c, F cmp) +{ + T y = as(c); + Function(y.begin(), y.end(), cmp); + f.col(i) = y; +} + +template ::type, + typename std::remove_reference::type>> + Function, + class FF, + class F> +void setResultParallelSection(mat &f, FF s, const int i, F cmp) +{ + T y; +#pragma omp critical + { + HELPER yy(s); + y = as(yy); + } + Function(y.begin(), y.end(), cmp); + f.col(i) = y; +} + +template Function, class FF> +void setResultParallelSection(mat &f, FF s, const int i) +{ + T y; +#pragma omp critical + { + HELPER yy(s); + y = as(yy); + } + Function(y.begin(), y.end()); + f.col(i) = y; +} + +template Function, class FF> +void setResult(mat &f, const int i, FF c) +{ + T y = as(c); + Function(y.begin(), y.end()); + f.col(i) = y; +} + +template ::type, + typename std::remove_reference::type>> + Function, + class FF, + class F> +void setResult(List &f, const int i, FF c, F cmp) +{ + T y = clone(as(c)); + Function(y.begin(), y.end(), cmp); + f[i] = y; +} + +template ::type, + typename std::remove_reference::type>> + Function, + class FF, + class F> +void setResultParallelSection(List &f, FF s, F cmp) +{ + T y; + int i; +#pragma omp critical + { + HELPER yy = *s; + y = as(yy); + i = s - f.begin(); + } + Function(y.begin(), y.end(), cmp); +#pragma omp critical + { + f[i] = HELPER(y.begin(), y.end()); + } +} + +template Function, class FF> +void setResultParallelSection(List &f, FF s) +{ + T y; + int i; +#pragma omp critical + { + HELPER yy = *s; + y = as(yy); + i = s - f.begin(); + } + Function(y.begin(), y.end()); +#pragma omp critical + { + f[i] = HELPER(y.begin(), y.end()); + } +} + +template Function, class FF> +void parallelSingleIteratorWithoutCopy(List &f, FF s) +{ + T y; + int i; +#pragma omp critical + { + HELPER yy(s); + y = as(yy); + i = s - f.begin(); + } + auto value = *Function(y.begin(), y.end()); +#pragma omp critical + { + f[i] = HELPER::create(value); + } +} + +template Function, class FF> +void parallelSingleIteratorWithoutCopy(colvec &f, FF s) +{ + T y; + int i; +#pragma omp critical + { + HELPER yy(s); + y = T(yy.begin(), yy.size(), false); + i = s - f.begin(); + } + auto value = *Function(y.begin(), y.end()); +#pragma omp critical + { + f[i] = value; + } +} + +template , typename T::iterator, typename T::iterator> Function, class FF> +RET parallelSingleIteratorWithoutCopy(FF s) +{ + T y; +#pragma omp critical + { + HELPER yy(s); + y = T(yy.begin(), yy.size(), false); + } + auto v = Function(y.begin(), y.end()); + return {static_cast(*v.first), static_cast(*v.second)}; +} + +template , typename T::iterator, typename T::iterator> Function, class FF> +RET singleIteratorWithoutCopy(FF c) +{ + Helper h(c); + T y(h.begin(), h.size(), false); + auto v = Function(y.begin(), y.end()); + return {static_cast(*v.first), static_cast(*v.second)}; +} + +inline long long int get_current_nanoseconds() +{ + return std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); +} +#endif \ No newline at end of file diff --git a/inst/include/Rfast2/matrix.hpp b/inst/include/Rfast2/matrix.hpp new file mode 100644 index 0000000..63505cc --- /dev/null +++ b/inst/include/Rfast2/matrix.hpp @@ -0,0 +1,171 @@ + +#ifndef MATRIX2_HPP +#define MATRIX2_HPP + +#include +#include +#include +#include +#include +#ifdef _OPENMP +#include +#endif + +#include "templates.h" +#include "helpers.hpp" +#include "vector.hpp" + +namespace Rfast +{ + using namespace Rcpp; + using namespace arma; + using namespace std; + using namespace chrono; + + inline NumericVector colTrimMean(NumericMatrix X, const double a = 0.05, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + mat x(X.begin(), X.nrow(), X.ncol(), false); + NumericVector f(x.n_cols); + colvec ff(f.begin(), f.size(), false); +#ifdef _OPENMP +#pragma omp parallel for if (parallel) num_threads(cores) +#endif + for (unsigned int i = 0; i < x.n_cols; ++i) + { + ff(i) = Rfast::TrimMean(x.col(i), a); + } + return f; + } + + inline NumericVector colTrimMean(DataFrame x, const double a = 0.05, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + NumericVector f(x.ncol()); + colvec ff(f.begin(), f.size(), false); + if (parallel) + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (DataFrame::iterator s = x.begin(); s < x.end(); ++s) + { + colvec y; + int i; +#ifdef _OPENMP +#pragma omp critical +#endif + { + NumericVector yy; + yy = *s; + y = colvec(yy.begin(), yy.size(), false); + i = s - x.begin(); + } + ff[i] = Rfast::TrimMean(y, a); + } + } + else + { + int i = 0; + NumericVector y(x.nrows()); + colvec yy; + for (auto c : x) + { + y = c; + yy = colvec(y.begin(), y.size(), false); + ff[i++] = Rfast::TrimMean(yy, a); + } + } + f.names() = as(x.names()); + return f; + } + + inline NumericVector rowTrimMean(NumericMatrix X, const double a = 0.05, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + mat x(X.begin(), X.nrow(), X.ncol(), false); + NumericVector f(x.n_rows); + colvec ff(f.begin(), f.size(), false); + +#ifdef _OPENMP +#pragma omp parallel for if (parallel) num_threads(cores) +#endif + for (unsigned int i = 0; i < x.n_rows; ++i) + { + ff(i) = TrimMean(x.row(i), a); + } + + return f; + } + + inline NumericMatrix colQuantile(DataFrame x, NumericVector Probs, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + colvec probs(Probs.begin(), Probs.size(), false); + NumericMatrix f(probs.n_elem, x.ncol()); + mat ff(f.begin(), probs.n_elem, x.ncol(), false); + if (parallel) + { +#ifdef _OPENMP +#pragma omp parallel for if (parallel) num_threads(cores) +#endif + for (DataFrame::iterator s = x.begin(); s < x.end(); ++s) + { + colvec y; + int i; +#ifdef _OPENMP +#pragma omp critical +#endif + { + NumericVector yy; + yy = *s; + y = colvec(yy.begin(), yy.size(), false); + i = s - x.begin(); + } + ff.col(i) = Quantile(y, probs); + } + } + else + { + int i = 0; + NumericVector y(x.nrows()); + colvec yy; + for (auto c : x) + { + y = c; + yy = colvec(y.begin(), y.size(), false); + ff.col(i++) = Quantile(yy, probs); + } + } + colnames(f) = as(x.names()); + return f; + } + + inline mat colQuantile(NumericMatrix X, NumericVector Probs, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + mat x(X.begin(), X.nrow(), X.ncol(), false); + colvec probs(Probs.begin(), Probs.size(), false); + mat f(probs.n_elem, x.n_cols); +#ifdef _OPENMP +#pragma omp parallel for if (parallel) num_threads(cores) +#endif + for (unsigned int i = 0; i < f.n_cols; ++i) + { + f.col(i) = Quantile(x.col(i), probs); + } + return f; + } + + inline mat rowQuantile(NumericMatrix X, NumericVector Probs, const bool parallel = false, const unsigned int cores = get_num_of_threads()) + { + mat x(X.begin(), X.nrow(), X.ncol(), false); + colvec probs(Probs.begin(), Probs.size(), false); + mat f(x.n_rows, probs.n_elem); +#ifdef _OPENMP +#pragma omp parallel for if (parallel) num_threads(cores) +#endif + for (unsigned int i = 0; i < f.n_rows; ++i) + { + f.row(i) = Quantile(x.row(i), probs); + } + return f; + } +} + +#endif diff --git a/inst/include/Rfast2/parallel.h b/inst/include/Rfast2/parallel.h new file mode 100644 index 0000000..c5bee7b --- /dev/null +++ b/inst/include/Rfast2/parallel.h @@ -0,0 +1,89 @@ + + +#ifndef PARALLEL_H +#define PARALLEL_H + +//[[Rcpp::depends(RcppParallel)]] +#include + +#include + +/* definition to expand macro then apply to pragma message */ + +// #define STRS(x) #x +// #define STR(x) STRS(x) +// #pragma message STR(__cplusplus) + +#if __cplusplus >= 201603L && (!defined(__APPLE__) && !defined(__MACH__) && !defined(__clang__)) +#define _PARALLEL_ +// #pragma message "Parallel is supported" +#include +#else +// #pragma message "Parallel is not supported" +#include +#endif + +namespace Rfast +{ + inline constexpr bool isStdParallelSupported() + { +#ifdef _PARALLEL_ + return true; +#else + return false; +#endif + } + + template + void sort(T begin, T end, const bool parallel = false) + { + if (parallel) + { +#ifdef _PARALLEL_ + std::sort(std::execution::par, begin, end); +#else + throw std::runtime_error("The C++ parallel library isn't supported by your system. Please, don't use the parallel argument."); +#endif + } + else + { + std::sort(begin, end); + } + } + + template + void nth_element(T begin, T middle, T end, const bool parallel = false) + { + if (parallel) + { +#ifdef _PARALLEL_ + std::nth_element(std::execution::par, begin, middle, end); +#else + throw std::runtime_error("The C++ parallel library isn't supported by your system. Please, don't use the parallel argument."); +#endif + } + else + { + std::nth_element(begin, middle, end); + } + } + + template + void nth_element(T begin, T middle, T end, Function cmp, const bool parallel = false) + { + if (parallel) + { +#ifdef _PARALLEL_ + std::nth_element(std::execution::par, begin, middle, end, cmp); +#else + throw std::runtime_error("The C++ parallel library isn't supported by your system. Please, don't use the parallel argument."); +#endif + } + else + { + std::nth_element(begin, middle, end, cmp); + } + } +} + +#endif \ No newline at end of file diff --git a/inst/include/Rfast2/templates.h b/inst/include/Rfast2/templates.h new file mode 100644 index 0000000..93e747a --- /dev/null +++ b/inst/include/Rfast2/templates.h @@ -0,0 +1,1294 @@ + +#ifndef TEMPLATES_H +#define TEMPLATES_H + +// [[Rcpp::depends(RcppArmadillo)]] +#include +#include +#include +#include +#include "parallel.h" + +using namespace std; +using namespace arma; +using namespace Rcpp; + +template +struct pr{ + f first; + s second; + bool is_good; + pr(f first=0,s second=0):first(first),second(second),is_good(false){} +}; + +typedef double (*Unary_Function)(double); // unary function +typedef double (*Binary_Function)(double,double); // binary function +typedef double (*Binary_Function_mat)(mat,double); // binary function + +template +using Mfunction = RET(*)(Args...); + +template +using ConditionFunction = bool(*)(T); + +// T: any simple data type. +template +inline bool notNA(T v){ + return !R_IsNA(v); +} + +/* +* F: unary function +* T1,T2: arguent class +* +* fill memory that starts from "startf" applying function F from "start" to "end" +*/ +template +void fill_with(T1 start,T1 end,T2 startf){ + for(;start!=end;++start,++startf){ + *startf=F(*start); + } +} + + +/* +* F: unary function +* T1,T2: arguent class +* +* fill memory that starts from "startf" applying function F to "x" +*/ +template +void fill_with(T1 x,T2 startf){ + for(typename T1::iterator start=x.begin();start!=x.end();++start){ + *startf=F(*start); + } +} + +/* +* F: unary function +* T1,T2: arguent class +* +* fill memory that starts from "startf" applying function F to "x" +*/ +template COND,typename T1,typename T2> +void fill_with_condition(T1 start,T1 end,T2 startf){ + for(;start!=end;++start,++startf){ + if(COND(*start)){ + *startf=F(*start); + } + } +} + + +/* +* T: argument and return class +* +* applying function F to "x" and return "x". +*/ +template +T foreach(T x){ + for(typename T::iterator start=x.begin();start!=x.end();++start){ + *start=F(*start); + } + return x; +} + + +/* +* T1,T2: argument class +* +* count the frequency of "value" in "x" +*/ +template +int count_value_helper(T1 x,T2 value){ + int s=0; + for(typename T1::iterator start=x.begin();start!=x.end();++start){ + if(*start==value){ + ++s; + } + } + return s; +} + + +/* +* T: argument class +* +* count the frequency of "NAs" in "x" +*/ +template +int count_NAs(T x){ + int s=0; + for(typename T::iterator start=x.begin();start!=x.end();++start){ + if(R_IsNA(*start)){ + ++s; + } + } + return s; +} + + +/* +* T: argument class +*/ +template +double med_helper(typename T::iterator start,typename T::iterator last){ + double F; + const int sz=last-start,middle=sz/2-1; + if(sz%2==0){ + std::nth_element(start,start+middle,last); + F=(start[middle]+*(std::min_element(start+middle+1,last)))/2.0; + }else{ + std::nth_element(start,start+middle+1,last); + F=start[middle+1]; + } + return F; +} + + +/* +* Ret: return class +* T: argument class +*/ +template +Ret Order(T x,const bool stable,const bool descend,const int init_v){ + Ret ind(x.size()); + iota(ind.begin(),ind.end(),init_v); + if(descend){ + auto descend_func = [&](int i,int j){return x[i-init_v]>x[j-init_v];}; + stable ? stable_sort(ind.begin(),ind.end(),descend_func) : sort(ind.begin(),ind.end(),descend_func); + }else{ + auto func = [&](int i,int j){return x[i-init_v] +void min_max(T *start,T *end,T &mn, T &mx){ + T xxx; + mn=mx=*start; + start++; + for(;start!=end;++start){ + xxx=*start; + if(xxx>mx){ + mx=xxx; + }else if(xxx +void maximum(T *start,T *end,T &mx){ + double xxx; + mx=*start; + start++; + for(;start!=end;++start){ + xxx=*start; + if(xxx>mx){ + mx=xxx; + } + } +} + + +/* +* T: argument class (not classes) +*/ +template +void minimum(T *start,T *end,T &mn){ + double xxx; + mn=*start; + start++; + for(;start!=end;++start){ + xxx=*start; + if(xxx +Ret Order_rank(T& x,const bool descend,const bool stable,const int n,const int k){ + Ret ind(x.size()-k); + iota(ind.begin(),ind.end(),0); + if(descend){ + auto descend_func = [&](int i,int j){return x[i]>x[j];}; + stable ? stable_sort(ind.begin(),ind.end()-n,descend_func) : sort(ind.begin(),ind.end()-n,descend_func); + }else{ + auto func = [&](int i,int j){return x[i] COND,class F> +T sum_with_condition(F x){ + T a=0; + for(typename F::iterator start=x.begin();start!=x.end();++start){ + if(COND(*start)){ + a+=*start; + } + } + return a; +} + +/* +* F: unary function +* T1,T2: arguent class +* +* fill memory that starts from "startf" applying function F from "start" to "end" +*/ +template COND,class F> +T sum_with_condition(F start,F end){ + T a=0; + for(;start!=end;++start){ + if(COND(*start)){ + a+=*start; + } + } + return a; +} + +/* +* F: a unary function +* T: argument class +* +* applying function F to "x" and sum the elements of "x" +*/ +template +double sum_with(T x){ + double a=0; + for(typename T::iterator start=x.begin();start!=x.end();++start){ + a+=F(*start); + } + return a; +} + +/* +* T: argument class +* +* applying function F to "x" and sum the elements of "x" +*/ +template F> +double sum_with(T *start,T *end){ + double a=0; + for(;start!=end;++start){ + a+=F(*start); + } + return a; +} + +/* +* F: a binary function +* T: argument class +* +* applying function F to "x" and sum the elements of "x" +*/ +template +double sum_with(T x,const double p){ + double a=0; + for(typename T::iterator start=x.begin();start!=x.end();++start){ + a+=F(*start,p); + } + return a; +} + +/* +* F: a binary function +* T: argument class +* +* applying function F to "x" and sum the elements of "x" +*/ +template +double sum_with(T1 x,T2& y){ + double a=0; + typename T1::iterator startx=x.begin(); + typename T2::iterator starty=y.begin(); + for(;startx!=x.end();++startx,++starty){ + a+=F(*startx,*starty); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x","y" and applying F2 to the result +*/ +template +double Apply(T1 x,T2& y,Binary_Function F1,Binary_Function F2){ + double a=0; + typename T1::iterator startx=x.begin(); + typename T2::iterator starty=y.begin(); + for(;startx!=x.end();++startx,++starty){ + a=F2(a,F1(*startx,*starty)); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x","y" and applying F2 to the result +*/ +template +double Apply(T1 x,T2& y){ + double a=0; + typename T1::iterator startx=x.begin(); + typename T2::iterator starty=y.begin(); + for(;startx!=x.end();++startx,++starty){ + a=F2(a,F1(*startx,*starty)); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x","y" and applying F2 to the result +*/ +template +double Apply(T1 *x,T1 *endx,T2 *y){ + double a=0; + for(;x!=endx;++x,++y){ + a=F2(a,F1(*x,*y)); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x","y" and applying F2 to the result +*/ +template +double Apply(T1 x,T2& y){ + double a=0; + typename T1::iterator startx=x.begin(); + typename T2::iterator starty=y.begin(); + for(;startx!=x.end();++startx,++starty){ + a=F3(a,F2(F1(*startx,*starty))); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x","y" and applying F2 to the result +*/ +template +double Apply(T1 x,T2& y,Binary_Function F1,Unary_Function F2,Binary_Function F3){ + double a=0; + typename T1::iterator startx=x.begin(); + typename T2::iterator starty=y.begin(); + for(;startx!=x.end();++startx,++starty){ + a=F3(a,F2(F1(*startx,*starty))); + } + return a; +} + +/* +* F1: a binary function +* F2: a binary function +* T1: argument class +* T2: argument class +* +* applying function F1 to "x" and applying F2 to the result +*/ +template +double Apply(T x){ + double a=0; + typename T::iterator start=x.begin(); + for(;start!=x.end();++start){ + a=F2(a,F1(*start)); + } + return a; +} + +/* +* T: return and argument type (Rcpp and simple types. Not armadillo) +*/ +template +T square2(T x){ + return x*x; +} + +/* +* Ret: return class +* T: argument class +*/ +template +Ret Tabulate(T x,int nroww){ + Ret f(nroww); + std::fill(f.begin(),f.end(),0); + typename Ret::iterator F=f.begin(); + typename T::iterator a=x.begin(); + for(;a!=x.end();++a){ + ++F[*a-1]; + } + return f; +} + +/* +* Ret: return class +* T: argument class +*/ +template +Ret Tabulate(int* start,int* end,int& nroww){ + Ret f(nroww); + std::fill(f.begin(),f.end(),0); + typename Ret::iterator F=f.begin(); + for(;start!=end;++start){ + ++F[*start-1]; + } + return f; +} + +/* +* Ret: return class +* T1: argument class +* T2: argument class +*/ +/* +* Ret: return class +* T1: argument class +* T2: argument class +*/ +template +Ret group_sum_helper(T1 x,T2 key,int *minn=nullptr,int *maxx=nullptr){ + int mn,mx; + const bool is_mn_null=(minn==nullptr),is_mx_null=(maxx==nullptr); + if(is_mx_null && is_mn_null){ + min_max(key.begin(),key.end(),mn,mx); + }else if(is_mx_null){ + mn=*minn; + maximum(key.begin(),key.end(),mx); + }else if(is_mn_null){ + mx=*maxx; + minimum(key.begin(),key.end(),mn); + }else{ + mn=*minn; + mx=*maxx; + } + typename T2::iterator kk=key.begin(); + vector f(mx-mn+1); + vector is_good(mx-mn+1); + typename T1::iterator xx=x.begin(),rr; + vector::iterator ff=f.begin(); + vector::iterator ok; + int index; + for(;xx!=x.end();++xx,++kk){ + index = *kk-mn; + is_good[index]=true; + f[index] += *xx; + } + int number_of_values=0; + for(auto v : is_good){ + if(v){ + ++number_of_values; + } + } + Ret res(number_of_values); + for(rr=res.begin(),ok=is_good.begin(),ff=f.begin();ff!=f.end();++ff){ + if(*ok++){ + *rr++=*ff; + } + } + return res; +} + + +/* Distixos gia tora to T prepei na einai panta mat +* dioti mono ta class tou arma exoun n_rows,n_cols +* Ret: return class +* T: argument class +*/ +template +Ret cross_x_y(T x,T y){ + const int ncl=x.n_cols,nrw=x.n_rows,p=y.n_cols; + Ret f(ncl,p); + Vec yv(nrw); + int i,j; + for(i=0;i +Ret cross_x(T x){ + const int ncl=x.n_cols; + Ret f(ncl,ncl); + double a; + int i,j; + for(i=0;i +Ret design_matrix_helper(T x) { + int i=0; + const int n=x.size(); + T tmp=sort_unique(x); + typename T::iterator xx=x.begin(),leksi_bg,leksi_en; + Ret Final(n,tmp.size()); + std::fill(Final.begin(),Final.end(),0); + for(leksi_bg=tmp.begin(),leksi_en=tmp.end(),i=0;xx!=x.end();++xx,++i){ + Final(i,lower_bound(leksi_bg,leksi_en,*xx)-leksi_bg)=1; + } + return Final; +} + +template +inline T madd(T x,T y){ + return x+y; +} + + +template +inline T mless(T x,T y){ + return x +inline T mless_eq(T x,T y){ + return x<=y; +} + +template +inline T mgreater_eq(T x,T y){ + return x>=y; +} + +template +inline T mgreater(T x,T y){ + return x>y; +} + + +template +inline T mdiv(T x,T y){ + return x/y; +} + + +template +inline T mdiff(T x,T y){ + return x-y; +} + + +template +inline T mmult(T x,T y){ + return x*y; +} + +template +inline T mequal(T x,T y){ + return x==y; +} + +template +inline T mmax(T x,T y){ + return std::max(x,y); +} + +template +inline T mmin(T x,T y){ + return std::min(x,y); +} + +template +inline T mand(T x,T y){ + return x&&y; +} + +template +inline T mor(T x,T y){ + return x||y; +} + +template +inline RET madd(T x,T y){ + return x+y; +} + + +template +inline RET mless(T x,T y){ + return x +inline RET mless_eq(T x,T y){ + return x<=y; +} + +template +inline RET mgreater_eq(T x,T y){ + return x>=y; +} + +template +inline RET mgreater(T x,T y){ + return x>y; +} + +template +inline bool mgreater(T x,T y){ + return x>y; +} + + +template +inline RET mdiv(T x,T y){ + return x/y; +} + + +template +inline RET mdiff(T x,T y){ + return x-y; +} + + +template +inline RET mmult(T x,T y){ + return x*y; +} + +template +inline RET mequal(T x,T y){ + return x==y; +} + +template +inline RET mmax(T x,T y){ + return std::max(x,y); +} + +template +inline RET mmin(T x,T y){ + return std::min(x,y); +} + +template +inline RET mand(T x,T y){ + return x&&y; +} + +template +inline RET mor(T x,T y){ + return x||y; +} + +template +inline RET madd(T1 x,T2 y){ + return x+y; +} + + +template +inline RET mless(T1 x,T2 y){ + return x +inline RET mless_eq(T1 x,T2 y){ + return x<=y; +} + +template +inline RET mgreater_eq(T1 x,T2 y){ + return x>=y; +} + +template +inline RET mgreater(T1 x,T2 y){ + return x>y; +} + + +template +inline RET mdiv(T1 x,T2 y){ + return x/y; +} + + +template +inline RET mdiff(T1 x,T2 y){ + return x-y; +} + + +template +inline RET mmult(T1 x,T2 y){ + return x*y; +} + +template +inline RET mequal(T1 x,T2 y){ + return x==y; +} + +template +inline RET mmax(T1 x,T2 y){ + return std::max(x,y); +} + +template +inline RET mmin(T1 x,T2 y){ + return std::min(x,y); +} + +template +inline RET mand(T1 x,T2 y){ + return x&&y; +} + +template +inline RET mor(T1 x,T2 y){ + return x||y; +} + + +template +void myoperator(T f[],T &x,T *y,int &len){ + T *ff=f; + for(int i=0;i +void as_integer_h_sorted(vector x,IntegerVector &f,const int init,const T val){ + const int n=x.size(); + int i,j=0,c=init; + sort(x.begin(),x.end()); + auto v=x[j]; + f[0]=init; + for(i=1;i +void as_integer_h(vector x,IntegerVector &f,const int init,const T val){ + const int n=x.size(); + int i,j=0,c=init; + vector ind=Order< vector,vector >(x,false,false,0); // diorthoseiii + x.push_back(val); + T v=x[ind[j]]; + f[ind[0]]=init; + for(i=1;i +void as_integer_h_with_names(vector x,List &L,const int init,const T val){ + const int n=x.size()+1; + int i,j=0,c=init; + vector ind=Order< vector,vector >(x,false,false,0); // diorthoseiii + x.push_back(val); + ind.push_back(0); + vector w; + T v=x[ind[j]]; + IntegerVector f(n-1); + f[ind[0]]=init; + for(i=1;i +vector table_use_na(vector x,const int use_na){ + typename vector::iterator new_end = std::remove_if(x.begin(),x.end(),R_IsNA); + sort(x.begin(),new_end); + vector n; + typename vector::iterator a=x.begin(),b=a+1; + List l; + int s=1; + for(;b!=new_end;++b){ + if(*a!=*b){ + n.push_back(s); + a=b; + s=1; + }else + ++s; + } + if(use_na==1){ + n.push_back(new_end-x.begin()); + } + return n; +} + +template +vector table_simple(vector x){ + sort(x.begin(),x.end()); + x.push_back(T(0)); + vector n; + typename vector::iterator a=x.begin(),b=a+1; + int s=1; + for(;b!=x.end();++b){ + if(*a!=*b){ + n.push_back(s); + a=b; + s=1; + }else{ + ++s; + } + } + return n; +} + +template +void table2_like_r(vector x,vector y,IntegerMatrix &f,const T val){ + const int n=x.size(); + int mx_x,mx_y; + IntegerVector ix(n),iy(n); + as_integer_h(x,ix,0,val); + as_integer_h(y,iy,0,val); + maximum(ix.begin(),ix.end(),mx_x); + maximum(iy.begin(),iy.end(),mx_y); + f=IntegerMatrix(mx_x+1,mx_y+1); + for(int i=0;i +void table2_like_r_with_names(vector x,vector y,List &L,const T val){ + const int n=x.size(); + int mx_x,mx_y; + List lx,ly; + as_integer_h_with_names(x,lx,0,val); + as_integer_h_with_names(y,ly,0,val); + IntegerVector ix=lx["f"],iy=ly["f"]; + maximum(ix.begin(),ix.end(),mx_x); + maximum(iy.begin(),iy.end(),mx_y); + IntegerMatrix f(mx_x+1,mx_y+1); + for(int i=0;i +Ret rank_mean(T x,const bool descend){ + const int n=x.size(),n_1=n+1; + int i,j=0; + x.resize(n_1); + x[n]=std::numeric_limits::max(); + I ind=Order_rank(x,descend,false,1,0); + int k=0,m,times=0; + Ret f(n); + double mn=0.0,v=x[ind[j]]; + for(i=1;i +Ret rank_max(T x,const bool descend){ + const int n=x.size(),n_1=n+1; + int i,j=0; + x.resize(n_1); + x[n]=std::numeric_limits::max(); + I ind=Order_rank(x,descend,false,1,0); + int k=0,m,times=0; + Ret f(n); + double v=x[ind[j]]; + for(i=1;i +Ret rank_min(T x,const bool descend){ + const int n=x.size(),n_1=n+1; + int i,j=0; + x.resize(n_1); + x[n]=std::numeric_limits::max(); + I ind=Order_rank(x,descend,false,1,0); + Ret f(n); + double v=x[ind[j]]; + f[ind[0]]=1; + for(i=1;i +Ret rank_first(T x,const bool descend,const bool stable){ + const int n=x.n_elem; + I ind=Order_rank(x,descend,stable,0,1); + Ret f(n); + for(int i=0;i +NumericVector eachcol_apply_helper(NumericMatrix& x,NumericVector& y,SEXP ind = R_NilValue, const bool parallel=false){ + const bool is_ind_null = Rf_isNull(ind); + const int n = is_ind_null ? x.ncol() : LENGTH(ind); + NumericVector f(n); + mat xx(x.begin(),x.nrow(),x.ncol(),false); + colvec yy(y.begin(),y.size(),false),ff(f.begin(),f.size(),false); + if(is_ind_null){ + if(parallel){ + #pragma omp parallel for + for(int i=0;i(xx.col(i),yy); + } + }else{ + for(int i=0;i(xx.col(i),yy); + } + } + }else{ + IntegerVector indd(ind); + icolvec iind(indd.begin(),indd.size(),false); + if(parallel){ + #pragma omp parallel for + for(int i=0;i(xx.col(iind[i]-1),yy); + } + }else{ + for(int i=0;i(xx.col(iind[i]-1),yy); + } + } + } + return f; +} + +template +double apply_eachrow_helper(SEXP x,SEXP y){ + int ncol=Rf_ncols(x),nrow=Rf_nrows(x); + SEXP mat=Rf_duplicate(x); + double *xx=REAL(mat),*end=xx+ncol*nrow,*yy=REAL(y),y3,*x3,s=0; + for(;xx!=end;++yy){ + y3=*yy; + for(x3=xx,xx+=nrow;x3!=xx;++x3){ + s=func(s,oper(*x3,y3)); + } + } + return s; +} + +template +SEXP eachrow_helper(SEXP x,SEXP y){ + int ncol=Rf_ncols(x),nrow=Rf_nrows(x); + SEXP mat=PROTECT(Rf_allocMatrix(type,nrow,ncol)); + T *xx=(T *) DATAPTR(x),*xend=xx+ncol*nrow,*yy=(T *) DATAPTR(y),yvalue,*x3; + RETURN_TYPE *m=(RETURN_TYPE*)DATAPTR(mat); + for(;xx!=xend;++yy){ + yvalue=*yy; + for(x3=xx,xx+=nrow;x3!=xx;++x3,++m){ + *m=oper(*x3,yvalue); + } + } + UNPROTECT(1); + return mat; +} + +template +double sum_x_op_y(SEXP x,SEXP y){ + const int n=LENGTH(x); + double s=0,*startx=REAL(x),*endx=startx+n,*starty=REAL(y); + for(;startx!=endx;++startx,++starty){ + s=func(s,oper(*startx,*starty)); + } + return s; +} + +template +double sum_x_op_x(SEXP x){ + const int n=LENGTH(x); + double s=0,*startx=REAL(x),*endx=startx+n; + for(;startx!=endx;++startx){ + s=func(s,oper(*startx,*startx)); + } + return s; +} + +//sort_unique,len_sort_unique,sort_int +template +void max_neg_pos(T *start,T *end,T &mx,T &mn,bool &pos,bool &neg){ + mn=mx=*start; + T x; + for(;start!=end;++start){ + x=*start; + if(x<0){ + neg=true; + if(mn>x){ + mn=x; + } + }else{ + pos=true; + if(mx Cond,Mfunction Oper> +int Apply_helper(int *start,int *end,int& val){ + int t=0; + for(;start!=end;++start){ + if(Cond(*start,val)){ + t=Oper(t,*start); + } + } + return t; +} + +template // less or greater +NumericVector negative_or_positive(NumericVector &x){ + double v,val=x[0]; + for(auto xx=x.begin()+1;xx!=x.end();++xx){ + v=*xx; + if(less_or_greater(v,0)){ + if(min_or_max(v,val)){ + val=v; + } + } + } + return NumericVector::create(val); +} + +template // less or greater +NumericVector negative_or_positive_min_max(NumericVector &x){ + double v,mn=x[0],mx=mn; + for(auto xx=x.begin()+1;xx!=x.end();++xx){ + v=*xx; + if(less_or_greater(v,0)){ + if(vmx){ + mx=v; + } + } + } + return NumericVector::create(mn,mx); +} + +template +double nth_simple(T& x,const int& elem,const bool parallel = false){ + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.end(),parallel); + return x[elem-1]; +} + +template +double nth_simple(T& x,const int& elem,const bool& descend,const bool parallel = false){ + descend ? + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.end(),[&](double a,double b){return a>b;},parallel) + : + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.end(),parallel); + + return x[elem-1]; +} + +template +double nth_na_rm(T& x,const int& elem,const bool& descend,const bool parallel = false){ + const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); + descend ? + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.begin()+new_end,[&](double a,double b){return a>b;},parallel) + : + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.begin()+new_end,parallel); + return x[elem-1]; +} + +template +double nth_helper(T& x,const int elem,const bool descend,const bool na_rm,const bool parallel = false){ + return na_rm ? nth_na_rm(x,elem,descend,parallel) : nth_simple(x,elem,descend,parallel); +} + + +template +int nth_index_simple(T& x,const int& elem,const bool& descend,const bool parallel = false){ + IntegerVector ind=seq(1,x.size()); + descend ? + Rfast::nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1]>x[j-1];},parallel) + : + Rfast::nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1] +int nth_index_na_rm(T& x,const int& elem,const bool& descend,const bool parallel = false){ + const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); + IntegerVector ind= seq(1,new_end); + descend ? + Rfast::nth_element(ind.begin(),ind.begin()+((elemx[j-1];},parallel) + : + Rfast::nth_element(ind.begin(),ind.begin()+((elem +int nth_helper_index(T& x,const int elem,const bool descend,const bool na_rm,const bool parallel = false){ + return na_rm ? nth_index_na_rm(x,elem,descend,parallel) : nth_index_simple(x,elem,descend,parallel); +} + +template +T nth_simple_n_elems(T& x,const int& elem,const bool& descend,const bool parallel = false){ + descend ? + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.end(),[&](double a,double b){return a>b;},parallel) + : + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.end(),parallel); + + return x(span(0,elem-1)); +} + +template +T nth_na_rm_n_elems(T& x,const int& elem,const bool& descend,const bool parallel = false){ + const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); + if(elemb;},parallel) + : + Rfast::nth_element(x.begin(),x.begin()+elem-1,x.begin()+new_end,parallel); + } + return x(span(0,elem-1)); +} + +template +T nth_helper_n_elems(T& x,const int elem,const bool descend,const bool na_rm,const bool parallel = false){ + return na_rm ? nth_na_rm_n_elems(x,elem,descend,parallel) : nth_simple_n_elems(x,elem,descend,parallel); +} + + +template +T nth_index_simple_n_elems(T& x,const int& elem,const bool& descend,const bool parallel = false){ + vec ind= linspace(1,x.size(),x.size()); + descend ? + Rfast::nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1]>x[j-1];},parallel) + : + Rfast::nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1] +T nth_index_na_rm_n_elems(T& x,const int& elem,const bool& descend,const bool parallel = false){ + const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); + vec ind= linspace(1,new_end,new_end); + descend ? + Rfast::nth_element(ind.begin(),ind.begin()+((elemx[j-1];},parallel) + : + Rfast::nth_element(ind.begin(),ind.begin()+((elem +T nth_helper_index_n_elems(T& x,const int elem,const bool descend,const bool na_rm,const bool parallel = false){ + return na_rm ? nth_index_na_rm_n_elems(x,elem,descend,parallel) : nth_index_simple_n_elems(x,elem,descend,parallel); +} + +#endif diff --git a/inst/include/Rfast2/templates_rfast2.h b/inst/include/Rfast2/templates_rfast2.h new file mode 100644 index 0000000..320f30b --- /dev/null +++ b/inst/include/Rfast2/templates_rfast2.h @@ -0,0 +1,67 @@ +#include + +using namespace std; +using namespace arma; +using namespace Rcpp; + +#ifndef TEMPLATES_RFAST2_H +#define TEMPLATES_RFAST2_H + +#include "templates.h" + +template func, const int init_val = 0> +SEXP group_col_h(SEXP x, SEXP gr, const int length_unique) +{ + const int ncl = Rf_ncols(x), nrw = Rf_nrows(x); + SEXP f = PROTECT(Rf_allocMatrix(TYPEOF(x), length_unique, ncl)); + int *ggr = INTEGER(gr); + T *ff = (T *)DATAPTR(f), *xx = (T *)DATAPTR(x); + for (int j = 0; j < length_unique * ncl; ++j) + { + ff[j] = init_val; + } + for (int j = 0; j < ncl; ++j) + { + const int col_index_f = j * length_unique, col_index_x = j * nrw; + for (int i = 0; i < nrw; ++i) + { + int ind_gr = ggr[i] - 1; + ff[ind_gr + col_index_f] = func(ff[ind_gr + col_index_f], xx[i + col_index_x]); + } + } + UNPROTECT(1); + return f; +} + +template +SEXP group_col_med_h(SEXP x, SEXP gr, const int length_unique) +{ + const int ncl = Rf_ncols(x), nrw = Rf_nrows(x); + SEXP f = PROTECT(Rf_allocMatrix(TYPEOF(x), length_unique, ncl)); + int *ggr = INTEGER(gr); + T *ff = (T *)DATAPTR(f), *xx = (T *)DATAPTR(x); + vector> eachcol_mat(length_unique, vector()); + for (int j = 0; j < length_unique * ncl; ++j) + { + ff[j] = 0; + } + for (int j = 0; j < ncl; ++j) + { + const int col_index_f = j * length_unique, col_index_x = j * nrw; + for (int i = 0; i < nrw; ++i) + { + int ind_gr = ggr[i] - 1; + eachcol_mat[ind_gr].push_back(xx[i + col_index_x]); + } + for (int i = 0; i < length_unique; ++i) + { + vector &tmp = eachcol_mat[i]; + ff[i + col_index_f] = med_helper>(tmp.begin(), tmp.end()); + tmp.clear(); + } + } + UNPROTECT(1); + return f; +} + +#endif \ No newline at end of file diff --git a/inst/include/Rfast2/vector.hpp b/inst/include/Rfast2/vector.hpp new file mode 100644 index 0000000..26d5e3d --- /dev/null +++ b/inst/include/Rfast2/vector.hpp @@ -0,0 +1,100 @@ + +#ifndef VECTOR2_HPP +#define VECTOR2_HPP + +#include +#include +#ifdef _OPENMP +#include +#endif + +#include "templates.h" +#include "assertions.hpp" + +#include "parallel.h" + +namespace Rfast +{ + + using namespace Rcpp; + using namespace arma; + using std::remove_if; + using std::string; + + template + Ret Quantile(T x, F &probs, const bool parallel = false) + { + Assertion::is_iterable::check_concept(); + Assertion::has_size::check_concept(); + Assertion::is_iterable::check_concept(); + Assertion::has_size::check_concept(); + + const unsigned int nprobs = probs.size(); + + Ret f(nprobs); + + 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()); + 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]; + f[i] = a + (h - hf) * (x[hf] - a); + } + } + else + { + for (unsigned int i = 0; i < nprobs; ++i) + { + double h = (x.size() - 1) * probs[i] + 1; + int hf = h; + double a, b; + if (probs[i] > 0.5) + { + a = nth_simple(x, hf - 1, false, parallel); + b = *min_element(x.begin() + hf, x.end()); + } + else + { + b = nth_simple(x, hf, false, parallel); + a = *max_element(x.begin(), x.begin() + hf); + } + f[i] = a + (h - hf) * (b - a); + } + } + + return f; + } + + template + double TrimMean(T x, const double a = 0.5, const bool parallel = false) + { + Assertion::is_iterable::check_concept(); + Assertion::has_size::check_concept(); + + const int n = x.size(); + int b1 = a * n; + int b11 = std::ceil(b1); + b1 = (b1 == b11) ? b11 + 1 : b11; + + const double a1 = nth_simple(x, b1, false, parallel); + const double a2 = nth_simple(x, n - b1 + 1, false, parallel); + double s = 0; + int p = 0; + for (auto xx : x) + { + if (xx >= a1 && xx <= a2) + { + s += xx; + ++p; + } + } + return s / p; + } + +} +#endif \ No newline at end of file diff --git a/man/Quantile.Rd b/man/Quantile.Rd index a5c78de..7c3b3f1 100644 --- a/man/Quantile.Rd +++ b/man/Quantile.Rd @@ -13,17 +13,17 @@ Sample quantiles and col/row wise quantiles. } \usage{ -colQuantile(x,probs,parallel=FALSE) -\method{colQuantile}{matrix}(x,probs,parallel=FALSE) -\method{colQuantile}{data.frame}(x,probs,parallel=FALSE) -rowQuantile(x,probs,parallel=FALSE) -Quantile(x,probs) +colQuantile(x,probs,parallel=FALSE,cores=0) +\method{colQuantile}{matrix}(x,probs,parallel=FALSE,cores=0) +\method{colQuantile}{data.frame}(x,probs,parallel=FALSE,cores=0) +rowQuantile(x,probs,parallel=FALSE,cores=0) +Quantile(x,probs,parallel=FALSE) } \arguments{ \item{x}{ Numeric vector whose sample quantiles are wanted. NA and NaN values are not allowed in numeric vectors. -For the col/row versions a numerical matrix. +For the col/row versions a numerical matrix or data.frame. } \item{probs}{ Numeric vector of probabilities with values in [0,1], not missing values. @@ -32,6 +32,10 @@ Values up to 2e-14 outside that range are accepted and automatically moved to th \item{parallel}{ Do you want to do it in parallel, for column - row major, in C++? TRUE or FALSE. } +\item{cores}{ +Number of cores to use for parallelism. Valid only when argument parallel is set to TRUE. +Default value is 0 and it means the maximum supported cores. +} } \details{ @@ -49,9 +53,6 @@ Manos Papadakis. R implementation and documentation: Manos Papadakis \email{papadakm95@gmail.com}. } -%\note{ -%% ~~further notes~~ -%} \seealso{ \code{ \link{trim.mean} } } diff --git a/man/Rfast2-package.Rd b/man/Rfast2-package.Rd index a7ca3f8..d24285e 100644 --- a/man/Rfast2-package.Rd +++ b/man/Rfast2-package.Rd @@ -18,8 +18,8 @@ Note 3: In general, make sure you give the correct input, in order to get the co \tabular{ll}{ Package: \tab Rfast2\cr Type: \tab Package\cr -Version: \tab 0.1.4 \cr -Date: \tab 2023-02-14\cr +Version: \tab 0.1.5.1 \cr +Date: \tab 2023-07-21\cr License: \tab GPL-2\cr } } diff --git a/man/boot.student2.Rd b/man/boot.student2.Rd index aa14eb7..2a5279b 100644 --- a/man/boot.student2.Rd +++ b/man/boot.student2.Rd @@ -28,7 +28,7 @@ The number of bootstrap samples to use. We bootstrap Student's (Gosset's) t-test statistic and not the Welch t-test statistic. For the latter case see the "boot.ttest2" function in Rfast. The difference is that Gosset's test statistic assumes equaility of the variances, which if violated leads to inlfated type I errors. Bootstrap calibration though takes care of this issue. -As for the bootstrap calibration, instead of sampling B times from each sample, we sample \eqn{sqrt{B}} from each +As for the bootstrap calibration, instead of sampling B times from each sample, we sample \eqn{\sqrt{B}} from each of them and then take all pairs. Each bootstrap sample is independent of each other, hence there is no violation of the theory (Chatzipantsiou et al., 2019). } @@ -64,8 +64,7 @@ R implementation and documentation: Michail Tsagris \email{mtsagris@uoc.gr}. \examples{ x <- rexp(40, 4) y <- rbeta(50, 2.5, 7.5) -system.time(t.test(x, y, var.equal = TRUE) ) -system.time( a <- boot.student2(x, y, 9999) ) -a +t.test(x, y, var.equal = TRUE) +boot.student2(x, y, 9999) } diff --git a/man/covequal.Rd b/man/covequal.Rd index f3555ca..cd0f4ba 100644 --- a/man/covequal.Rd +++ b/man/covequal.Rd @@ -1,12 +1,12 @@ -\name{Hypothesis tests for equality of a covariance matrix} +\name{Hypothesis test for equality of a covariance matrix} \alias{covequal} \title{ -Hypothesis tests for equality of a covariance matrix +Hypothesis test for equality of a covariance matrix } \description{ -Hypothesis tests for equality of a covariance matrix. +Hypothesis test for equality of a covariance matrix. } \usage{ diff --git a/man/kumar.mle.Rd b/man/kumar.mle.Rd index 2fb8880..fb152ad 100644 --- a/man/kumar.mle.Rd +++ b/man/kumar.mle.Rd @@ -1,7 +1,7 @@ \name{MLE of distributions defined for proportions} \alias{kumar.mle} -\alias{zil.mle} \alias{simplex.mle} +\alias{zil.mle} \alias{unitweibull.mle} \alias{cbern.mle} \alias{sp.mle} @@ -49,9 +49,9 @@ Usually a list with three elements, but this is not for all cases. \item{iters}{The number of iterations required for the Newton-Raphson to converge. } \item{param}{ -The two parameters (shape and scale) of the Kumaraswamy distribution or the means -and sigma of the simpled distribution. For the zero inflated logistic normal, -the probability of non zeros, the mean and the unbiased variance. +The two parameters (shape and scale) of the Kumaraswamy distribution. +For the zero inflated logistic normal, the probability of non zeros, +the mean and the unbiased variance. } \item{loglik}{The value of the maximised log-likelihood. } @@ -64,11 +64,7 @@ random processes. Journal of Hydrology 46(1-2): 79-88. Jones M.C. (2009). Kumaraswamy's distribution: A beta-type distribution with some tractability advantages. Statistical Methodology, 6(1): 70-81. -Zhang W. & Wei H. (2008). Maximum likelihood estimation for simplex distribution nonlinear -mixed models via the stochastic approximation algorithm. -The Rocky Mountain Journal of Mathematics, 38(5): 1863-1875. - -Mazucheli J., Menezes A.F.B., Fernandes L.B., de Oliveira R.P. & Ghitany M.E. (2020). +Mazucheli J., Menezes A.F.B., Fernandes L.B., de Oliveira R.P. and Ghitany M.E. (2020). The unit-Weibull distribution as an alternative to the Kumaraswamy distribution for the modeling of quantiles conditional on covariates. Journal of Applied Statistics, DOI:10.1080/02664763.2019.1657813 diff --git a/man/stud.ttests.Rd b/man/stud.ttests.Rd index 16413be..d69db81 100644 --- a/man/stud.ttests.Rd +++ b/man/stud.ttests.Rd @@ -67,10 +67,10 @@ R implementation and documentation: Michail Tsagris \email{mtsagris@uoc.gr}. x = matrix( rnorm(100 * 20), ncol = 20) ## 100 observations in total ina = rbinom(100, 1, 0.6) + 1 ## independent samples t-test -system.time( stud.ttests(x, ina = ina) ) +stud.ttests(x, ina = ina) x1 = x[ina == 1, ] x2 = x[ina == 2, ] -system.time( stud.ttests(x1, x2) ) +stud.ttests(x1, x2) x <- NULL } diff --git a/man/trim.mean.Rd b/man/trim.mean.Rd index ee9e4f8..925676b 100644 --- a/man/trim.mean.Rd +++ b/man/trim.mean.Rd @@ -12,16 +12,16 @@ Trimmed mean. } \usage{ -trim.mean(x, a = 0.05) -colTrimMean(x, a = 0.05,parallel=FALSE) -\method{colTrimMean}{matrix}(x,a = 0.05,parallel=FALSE) -\method{colTrimMean}{data.frame}(x,a = 0.05,parallel=FALSE) -rowTrimMean(x, a = 0.05,parallel=FALSE) +trim.mean(x, a = 0.05,parallel=FALSE) +colTrimMean(x, a = 0.05,parallel=FALSE,cores=0) +\method{colTrimMean}{matrix}(x,a = 0.05,parallel=FALSE,cores=0) +\method{colTrimMean}{data.frame}(x,a = 0.05,parallel=FALSE,cores=0) +rowTrimMean(x, a = 0.05,parallel=FALSE,cores=0) } \arguments{ \item{x}{ -A numerical vector or a numerical matrix. +A numerical vector or a numerical matrix or data.frame. } \item{a}{ A number in (0, 1), the proportion of data to trim. @@ -29,6 +29,10 @@ A number in (0, 1), the proportion of data to trim. \item{parallel}{ Run the algorithm parallel in C++. } +\item{cores}{ +Number of cores to use for parallelism. Valid only when argument parallel is set to TRUE. +Default value is 0 and it means the maximum supported cores. +} } \details{ @@ -51,10 +55,6 @@ Michail Tsagris and Manos Papadakis. R implementation and documentation: Michail Tsagris \email{mtsagris@uoc.gr} and Manos Papadakis \email{papadakm95@gmail.com}. } -%\note{ -%% ~~further notes~~ -%} - \seealso{ \code{ \link{Quantile} } diff --git a/src/Eval.cpp b/src/Eval.cpp index bc042cd..b9e4e4d 100644 --- a/src/Eval.cpp +++ b/src/Eval.cpp @@ -4,7 +4,7 @@ #include #include #include -#include "templates.h" +#include "Rfast2/templates.h" #include using namespace std; diff --git a/src/Makevars b/src/Makevars index f9bb09c..d864e97 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,4 +1,4 @@ -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) `"${R_HOME}/bin/Rscript" -e "RcppParallel::LdFlags()"` PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/Makevars.win b/src/Makevars.win index 57a5534..e7695f2 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,4 +1,4 @@ -PKG_LIBS += $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) +PKG_LIBS += $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) `"${R_HOME}/bin/Rscript" -e "RcppParallel::LdFlags()"` PKG_CXXFLAGS= $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/Random.h b/src/Random.h index ebfad9c..5753af7 100644 --- a/src/Random.h +++ b/src/Random.h @@ -44,7 +44,7 @@ namespace Random { return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); } }; - }; + } template class uniform : public internal::Integer_Core{ @@ -106,6 +106,6 @@ namespace Random { } }; -}; +} #endif \ No newline at end of file diff --git a/src/apply_funcs_templates.h b/src/apply_funcs_templates.h index 9bdd3ae..affc673 100644 --- a/src/apply_funcs_templates.h +++ b/src/apply_funcs_templates.h @@ -1,4 +1,4 @@ -#include "templates.h" +#include "Rfast2/templates.h" using namespace std; diff --git a/src/censweib_reg.cpp b/src/censweib_reg.cpp index 9a2551c..709d7d7 100644 --- a/src/censweib_reg.cpp +++ b/src/censweib_reg.cpp @@ -2,7 +2,7 @@ #include #include "mn.h" -#include "templates.h" +#include "Rfast2/templates.h" // [[Rcpp::depends(RcppArmadillo)]] using namespace arma; diff --git a/src/col_row_utilities.cpp b/src/col_row_utilities.cpp index 7a8b6f9..25647a2 100644 --- a/src/col_row_utilities.cpp +++ b/src/col_row_utilities.cpp @@ -1,5 +1,6 @@ #include -#include "templates.h" +#include "Rfast2.h" +#include "Rfast2/templates_rfast2.h" using namespace Rcpp; @@ -29,248 +30,136 @@ void group_col_vars_h(SEXP& x,SEXP& gr,const int length_unique,Environment& resu UNPROTECT(2); }*/ -SEXP group_col(SEXP x,SEXP y,const int length_unique,const string method="sum"){ - if(method == "sum"){ - if(Rf_isInteger(x)) - return group_col_h>(x,y,length_unique); - else if(Rf_isReal(x)) - return group_col_h>(x,y,length_unique); +SEXP group_col(SEXP x, SEXP y, const int length_unique, const string method = "sum") +{ + if (method == "sum") + { + if (Rf_isInteger(x)) + return group_col_h>(x, y, length_unique); + else if (Rf_isReal(x)) + return group_col_h>(x, y, length_unique); else stop("Error: Unsupported type of matrix."); - }else if(method == "max"){ - if(Rf_isInteger(x)) - return group_col_h,INT_MIN>(x,y,length_unique); - else if(Rf_isReal(x)) - return group_col_h,INT_MIN>(x,y,length_unique); + } + else if (method == "max") + { + if (Rf_isInteger(x)) + return group_col_h, INT_MIN>(x, y, length_unique); + else if (Rf_isReal(x)) + return group_col_h, INT_MIN>(x, y, length_unique); else stop("Error: Unsupported type of matrix."); - }else if(method == "min"){ - if(Rf_isInteger(x)) - return group_col_h,INT_MAX>(x,y,length_unique); - else if(Rf_isReal(x)) - return group_col_h,INT_MAX>(x,y,length_unique); + } + else if (method == "min") + { + if (Rf_isInteger(x)) + return group_col_h, INT_MAX>(x, y, length_unique); + else if (Rf_isReal(x)) + return group_col_h, INT_MAX>(x, y, length_unique); + else + stop("Error: Unsupported type of matrix."); + } /*else if(method == "var"){ + if(Rf_isInteger(x)) + group_col_vars_h(x,y,length_unique,result); + else if(Rf_isReal(x)) + group_col_vars_h(x,y,length_unique,result); + else + stop("Error: Unsupported type of matrix."); + return R_NilValue; + }*/ + else if (method == "median") + { + if (Rf_isInteger(x)) + return group_col_med_h(x, y, length_unique); + else if (Rf_isReal(x)) + return group_col_med_h(x, y, length_unique); else stop("Error: Unsupported type of matrix."); - }/*else if(method == "var"){ - if(Rf_isInteger(x)) - group_col_vars_h(x,y,length_unique,result); - else if(Rf_isReal(x)) - group_col_vars_h(x,y,length_unique,result); - else - stop("Error: Unsupported type of matrix."); - return R_NilValue; - }*/else if(method == "median"){ - if(Rf_isInteger(x)) - return group_col_med_h(x,y,length_unique); - else if(Rf_isReal(x)) - return group_col_med_h(x,y,length_unique); - else - stop("Error: Unsupported type of matrix."); } stop("Error: Unsupported method.\n"); return R_NilValue; } - -RcppExport SEXP Rfast2_col_group(SEXP x,SEXP y,SEXP length_uniqueSEXP,SEXP methodSEXP){ -BEGIN_RCPP - RObject __result; - RNGScope __rngScope; - traits::input_parameter< const int >::type length_unique(length_uniqueSEXP); - traits::input_parameter< const string >::type method(methodSEXP); - __result = group_col(x,y,length_unique,method); - return __result; -END_RCPP -} - -/***********************************************************************************/ - -mat colQuantile(NumericMatrix X,NumericVector Probs,const bool parallel = false){ - mat x(X.begin(),X.nrow(),X.ncol(),false); - colvec probs(Probs.begin(),Probs.size(),false); - mat f(probs.n_elem,x.n_cols); - if(parallel){ - #pragma omp parallel for - for(unsigned int i=0;i(x.col(i),probs); - } - }else{ - for(unsigned int i=0;i(x.col(i),probs); - } - } - return f; +RcppExport SEXP Rfast2_col_group(SEXP x, SEXP y, SEXP length_uniqueSEXP, SEXP methodSEXP) +{ + BEGIN_RCPP + RObject __result; + RNGScope __rngScope; + traits::input_parameter::type length_unique(length_uniqueSEXP); + traits::input_parameter::type method(methodSEXP); + __result = group_col(x, y, length_unique, method); + return __result; + END_RCPP } -//[[Rcpp::export]] -NumericMatrix colQuantile(DataFrame x,NumericVector Probs,const bool parallel = false){ - colvec probs(Probs.begin(),Probs.size(),false); - NumericMatrix f(probs.n_elem,x.ncol()); - mat ff(f.begin(),probs.n_elem,x.ncol(),false); - if(parallel){ - #pragma omp parallel for - for(DataFrame::iterator s = x.begin(); s < x.end(); ++s){ - colvec y; - int i; - #pragma omp critical - { - NumericVector yy; - yy=*s; - y = colvec(yy.begin(),yy.size(),false); - i = s-x.begin(); - } - ff.col(i)=Quantile(y,probs); - } - }else{ - int i=0; - NumericVector y(x.nrows()); - colvec yy; - for(auto c : x){ - y=c; - yy = colvec(y.begin(),y.size(),false); - ff.col(i++)=Quantile(yy,probs); - } +RcppExport SEXP Rfast2_colQuantile(SEXP xSEXP, SEXP ProbsSEXP, SEXP parallelSEXP, SEXP coresSEXP) +{ + BEGIN_RCPP + RObject __result; + RNGScope __rngScope; + traits::input_parameter::type Probs(ProbsSEXP); + traits::input_parameter::type parallel(parallelSEXP); + traits::input_parameter::type cores(coresSEXP); + if (Rf_isNewList(xSEXP)) + { + DataFrame x(xSEXP); + __result = Rfast::colQuantile(x, Probs, parallel, cores); } - colnames(f) = as(x.names()); - return f; -} - -RcppExport SEXP Rfast2_colQuantile(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP){ -BEGIN_RCPP - RObject __result; - RNGScope __rngScope; - traits::input_parameter< NumericVector >::type Probs(ProbsSEXP); - traits::input_parameter< const bool >::type parallel(parallelSEXP); - if(Rf_isNewList(xSEXP)){ - DataFrame x(xSEXP); - __result = colQuantile(x,Probs,parallel); - }else if(Rf_isMatrix(xSEXP)){ - NumericMatrix x(xSEXP); - __result = colQuantile(x,Probs,parallel); - } - return __result; -END_RCPP -} - -/***********************************************************************************/ - -mat rowQuantile(NumericMatrix X,NumericVector Probs,const bool parallel){ - mat x(X.begin(),X.nrow(),X.ncol(),false); - colvec probs(Probs.begin(),Probs.size(),false); - mat f(x.n_rows,probs.n_elem); - if(parallel){ - #pragma omp parallel for - for(unsigned int i=0;i(x.row(i),probs); - } - }else{ - for(unsigned int i=0;i(x.row(i),probs); - } - } - return f; -} - -RcppExport SEXP Rfast2_rowQuantile(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP){ -BEGIN_RCPP - RObject __result; - RNGScope __rngScope; - traits::input_parameter< NumericMatrix >::type x(xSEXP); - traits::input_parameter< NumericVector >::type Probs(ProbsSEXP); - traits::input_parameter< const bool >::type parallel(parallelSEXP); - __result = rowQuantile(x,Probs,parallel); - return __result; -END_RCPP -} - -/************************************************************************************/ - -NumericVector colTrimMean(NumericMatrix X,const double a=0.05,const bool parallel=false){ - mat x(X.begin(),X.nrow(),X.ncol(),false); - NumericVector f(x.n_cols); - colvec ff(f.begin(),f.size(),false); - if(parallel){ -#pragma omp parallel for - for(unsigned int i=0;i(x.col(i),a); - }else{ - for(unsigned int i=0;i(x.col(i),a); - } - return f; -} - -NumericVector colTrimMean(DataFrame x,const double a=0.05,const bool parallel=false){ - NumericVector f(x.ncol()); - colvec ff(f.begin(),f.size(),false); - if(parallel){ - #pragma omp parallel for - for(DataFrame::iterator s = x.begin(); s < x.end(); ++s){ - colvec y; - int i; - #pragma omp critical - { - NumericVector yy; - yy=*s; - y = colvec(yy.begin(),yy.size(),false); - i = s-x.begin(); - } - ff[i]=trimmean_h(y,a); - } - }else{ - int i=0; - NumericVector y(x.nrows()); - colvec yy; - for(auto c : x){ - y=c; - yy = colvec(y.begin(),y.size(),false); - ff[i++]=trimmean_h(yy,a); - } + else if (Rf_isMatrix(xSEXP)) + { + NumericMatrix x(xSEXP); + __result = Rfast::colQuantile(x, Probs, parallel, cores); } - f.names() = as(x.names()); - return f; + return __result; + END_RCPP } -RcppExport SEXP Rfast2_colTrimMean(SEXP xSEXP,SEXP aSEXP,SEXP parallelSEXP){ -BEGIN_RCPP - RObject __result; - RNGScope __rngScope; - traits::input_parameter< const double >::type a(aSEXP); - traits::input_parameter< const bool >::type parallel(parallelSEXP); - if(Rf_isNewList(xSEXP)){ - DataFrame x(xSEXP); - __result = colTrimMean(x,a,parallel); - }else if(Rf_isMatrix(xSEXP)){ - NumericMatrix x(xSEXP); - __result = colTrimMean(x,a,parallel); - } - return __result; -END_RCPP +RcppExport SEXP Rfast2_rowQuantile(SEXP xSEXP, SEXP ProbsSEXP, SEXP parallelSEXP, SEXP coresSEXP) +{ + BEGIN_RCPP + RObject __result; + RNGScope __rngScope; + traits::input_parameter::type x(xSEXP); + traits::input_parameter::type Probs(ProbsSEXP); + traits::input_parameter::type parallel(parallelSEXP); + traits::input_parameter::type cores(coresSEXP); + __result = Rfast::rowQuantile(x, Probs, parallel, cores); + return __result; + END_RCPP } -NumericVector rowTrimMean(NumericMatrix X,const double a=0.05,const bool parallel=false){ - mat x(X.begin(),X.nrow(),X.ncol(),false); - NumericVector f(x.n_rows); - colvec ff(f.begin(),f.size(),false); - if(parallel){ -#pragma omp parallel for - for(unsigned int i=0;i(x.row(i),a); - }else{ - for(unsigned int i=0;i(x.row(i),a); - } - return f; +RcppExport SEXP Rfast2_colTrimMean(SEXP xSEXP, SEXP aSEXP, SEXP parallelSEXP, SEXP coresSEXP) +{ + BEGIN_RCPP + RObject __result; + RNGScope __rngScope; + traits::input_parameter::type a(aSEXP); + traits::input_parameter::type parallel(parallelSEXP); + traits::input_parameter::type cores(coresSEXP); + if (Rf_isNewList(xSEXP)) + { + DataFrame x(xSEXP); + __result = Rfast::colTrimMean(x, a, parallel, cores); + } + else if (Rf_isMatrix(xSEXP)) + { + NumericMatrix x(xSEXP); + __result = Rfast::colTrimMean(x, a, parallel, cores); + } + return __result; + END_RCPP } -RcppExport SEXP Rfast2_rowTrimMean(SEXP xSEXP,SEXP aSEXP,SEXP parallelSEXP){ -BEGIN_RCPP - RObject __result; - RNGScope __rngScope; - traits::input_parameter< NumericMatrix >::type X(xSEXP); - traits::input_parameter< const double >::type a(aSEXP); - traits::input_parameter< const bool >::type parallel(parallelSEXP); - __result = rowTrimMean(X,a,parallel); - return __result; -END_RCPP +RcppExport SEXP Rfast2_rowTrimMean(SEXP xSEXP, SEXP aSEXP, SEXP parallelSEXP, SEXP coresSEXP) +{ + BEGIN_RCPP + RObject __result; + RNGScope __rngScope; + traits::input_parameter::type X(xSEXP); + traits::input_parameter::type a(aSEXP); + traits::input_parameter::type parallel(parallelSEXP); + traits::input_parameter::type cores(coresSEXP); + __result = Rfast::rowTrimMean(X, a, parallel, cores); + return __result; + END_RCPP } diff --git a/src/colbeta_mle.cpp b/src/colbeta_mle.cpp index cf844c9..ae49570 100644 --- a/src/colbeta_mle.cpp +++ b/src/colbeta_mle.cpp @@ -2,7 +2,7 @@ #include #include "reg_lib2.h" -#include "templates.h" +#include "Rfast2/templates.h" #include "mn.h" // [[Rcpp::depends(RcppArmadillo)]] diff --git a/src/colcauchy_mle.cpp b/src/colcauchy_mle.cpp index e30fa5e..d149271 100644 --- a/src/colcauchy_mle.cpp +++ b/src/colcauchy_mle.cpp @@ -2,7 +2,7 @@ #include #include "reg_lib2.h" -#include "templates.h" +#include "Rfast2/templates.h" #include "mn.h" // [[Rcpp::depends(RcppArmadillo)]] diff --git a/src/colspml_mle.cpp b/src/colspml_mle.cpp index 9a3be66..f4f3c3d 100644 --- a/src/colspml_mle.cpp +++ b/src/colspml_mle.cpp @@ -2,7 +2,7 @@ #include #include -#include "templates.h" +#include "Rfast2/templates.h" #include "mn.h" #include "reg_lib2.h" diff --git a/src/cts.h b/src/cts.h index 564fc75..8bcff87 100644 --- a/src/cts.h +++ b/src/cts.h @@ -7,7 +7,7 @@ #include #include #include -#include "templates.h" +#include "Rfast2/templates.h" // [[Rcpp::plugins("cpp11")]] // [[Rcpp::depends("RcppArmadillo")]] diff --git a/src/gamma_reg.cpp b/src/gamma_reg.cpp index 5bbec9f..8096434 100644 --- a/src/gamma_reg.cpp +++ b/src/gamma_reg.cpp @@ -2,7 +2,7 @@ #include #include "reg_lib2.h" -#include "templates.h" +#include "Rfast2/templates.h" #include "apply_funcs_templates.h" // [[Rcpp::depends(RcppArmadillo)]] diff --git a/src/gamma_regs.cpp b/src/gamma_regs.cpp index 5f619a1..2894da2 100644 --- a/src/gamma_regs.cpp +++ b/src/gamma_regs.cpp @@ -2,7 +2,7 @@ #include #include "reg_lib2.h" -#include "templates.h" +#include "Rfast2/templates.h" #include "apply_funcs_templates.h" // [[Rcpp::depends(RcppArmadillo)]] diff --git a/src/init.c b/src/init.c index cd89c92..e579ca4 100644 --- a/src/init.c +++ b/src/init.c @@ -1,102 +1,99 @@ #include #include -//Manos +// Manos -SEXP Rfast2_benchmark(SEXP exprsSEXP,SEXP envSEXP,SEXP timSEXP,SEXP indicesSEXP); -SEXP Rfast2_col_group(SEXP x,SEXP y,SEXP length_uniqueSEXP,SEXP methodSEXP); -SEXP Rfast2_colQuantile(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP); -SEXP Rfast2_is_upper_tri(SEXP xSEXP,SEXP dgSEXP); -SEXP Rfast2_is_lower_tri(SEXP xSEXP,SEXP dgSEXP); -SEXP Rfast2_is_skew_symmetric(SEXP xSEXP); -SEXP Rfast2_lud(SEXP xSEXP); -SEXP Rfast2_merge(SEXP xSEXP,SEXP ySEXP); -SEXP Rfast2_mmpc2(SEXP ySEXP,SEXP xSEXP,SEXP max_kSEXP,SEXP thresholdSEXP,SEXP testSEXP,SEXP Ini,SEXP parallelSEXP,SEXP maxitersSEXP,SEXP tolSEXP,SEXP backwardSEXP); -SEXP Rfast2_Quantile(SEXP xSEXP,SEXP ProbsSEXP); -SEXP Rfast2_rowQuantile(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP); -SEXP Rfast2_colTrimMean(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP); -SEXP Rfast2_rowTrimMean(SEXP xSEXP,SEXP ProbsSEXP,SEXP parallelSEXP); -SEXP Rfast2_Runif(SEXP nSEXP,SEXP minSEXP,SEXP maxSEXP); -SEXP Rfast2_Sample_int(SEXP nSEXP,SEXP sizeSEXP,SEXP replaceSEXP); -SEXP Rfast2_Sample(SEXP xSEXP,SEXP sizeSEXP,SEXP replaceSEXP); -SEXP Rfast2_trimmean(SEXP xSEXP,SEXP aSEXP); +SEXP Rfast2_benchmark(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_col_group(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_colQuantile(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_is_upper_tri(SEXP, SEXP); +SEXP Rfast2_is_lower_tri(SEXP, SEXP); +SEXP Rfast2_is_skew_symmetric(SEXP); +SEXP Rfast2_lud(SEXP); +SEXP Rfast2_merge(SEXP, SEXP); +SEXP Rfast2_mmpc2(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_Quantile(SEXP, SEXP, SEXP); +SEXP Rfast2_rowQuantile(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_colTrimMean(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_rowTrimMean(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_Runif(SEXP, SEXP, SEXP); +SEXP Rfast2_Sample_int(SEXP, SEXP, SEXP); +SEXP Rfast2_Sample(SEXP, SEXP, SEXP); +SEXP Rfast2_trimmean(SEXP, SEXP, SEXP); -//Manos +// Manos -//Marios +// Marios -SEXP Rfast2_inter(SEXP xSEXP,SEXP ySEXP); -SEXP Rfast2_mmp_c(SEXP target_varsSEXP,SEXP dsSEXP,SEXP max_kSEXP,SEXP thresSEXP,SEXP methodSEXP,SEXP initsSEXP,SEXP hash_onSEXP,SEXP stats_kvSEXP,SEXP pvalues_kvSEXP,SEXP bws_onSEXP); +SEXP Rfast2_inter(SEXP, SEXP); +SEXP Rfast2_mmp_c(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -//Marios +// Marios -//Stefanos +// Stefanos -SEXP Rfast2_add_term(SEXP YSEXP, SEXP XincSEXP, SEXP XoutSEXP, SEXP devi_0SEXP,SEXP typeSEXP,SEXP tolSEXP,SEXP logged,SEXP parallel,SEXP maxiters); -SEXP Rfast2_colspml_mle(SEXP xSEXP, SEXP tolSEXP,SEXP maxitersSEXP,SEXP parallelSEXP); -SEXP Rfast2_colcauchy_mle(SEXP XSEXP,SEXP tolSEXP,SEXP parallelSEXP,SEXP maxitersSEXP); -SEXP Rfast2_colbeta_mle(SEXP XSEXP,SEXP tolSEXP,SEXP parallelSEXP,SEXP maxitersSEXP); -SEXP Rfast2_censweib_reg(SEXP YSEXP,SEXP XSEXP,SEXP diSEXP,SEXP tolSEXP,SEXP maxitersSEXP); -//SEXP Rfast2_frechet2_c(SEXP XSEXP, SEXP DiSEXP, SEXP aSEXP, SEXP k1SEXP); -SEXP Rfast2_fbed_reg(SEXP ySEXP,SEXP xSEXP,SEXP sigSEXP,SEXP typeSEXP,SEXP idSEXP,SEXP kSEXP,SEXP backwardSEXP, SEXP tolSEXP,SEXP parallelSEXP,SEXP maxitersSEXP); -SEXP Rfast2_fedhc_skeleton(SEXP XSEXP, SEXP INI_PVALSEXP, SEXP nSEXP,SEXP laSEXP,SEXP methodSEXP,SEXP RmatSEXP, SEXP parallelSEXP); -SEXP Rfast2_multinom_reg(SEXP ySEXP,SEXP x0SEXP, SEXP tolSEXP,SEXP maxitersSEXP); -SEXP Rfast2_weib_regs(SEXP ySEXP,SEXP xSEXP, SEXP tolSEXP,SEXP loggedSEXP,SEXP maxitersSEXP, SEXP parallelSEXP); -SEXP Rfast2_welch_tests(SEXP xSEXP, SEXP ySEXP, SEXP loggedSEXP, SEXP parallelSEXP); -SEXP Rfast2_wild_boot(SEXP xSEXP, SEXP ySEXP, SEXP clusterSEXP, SEXP indSEXP, SEXP RSEXP, SEXP tabSEXP, SEXP parallelSEXP); -SEXP Rfast2_mmhc_skeleton(SEXP XSEXP, SEXP INI_PVALSEXP, SEXP nSEXP,SEXP laSEXP,SEXP maxkSEXP,SEXP methodSEXP,SEXP RmatSEXP, SEXP parallelSEXP); -SEXP Rfast2_negbin_reg(SEXP ySEXP, SEXP xSEXP, SEXP tolSEXP, SEXP maxitersSEXP); -SEXP Rfast2_negbin_regs(SEXP ySEXP, SEXP xSEXP, SEXP tolSEXP, SEXP maxitersSEXP, SEXP parallelSEXP); -SEXP Rfast2_gamma_regs(SEXP YSEXP, SEXP XSEXP,SEXP tolSEXP,SEXP loggedSEXP,SEXP parallelSEXP,SEXP maxitersSEXP); -SEXP Rfast2_gamma_reg(SEXP YSEXP, SEXP XSEXP,SEXP modSEXP,SEXP tolSEXP,SEXP maxitersSEXP); +SEXP Rfast2_add_term(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_colspml_mle(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_colcauchy_mle(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_colbeta_mle(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_censweib_reg(SEXP, SEXP, SEXP, SEXP, SEXP); +// SEXP Rfast2_frechet2_c(SEXP , SEXP , SEXP , SEXP k1SEXP); +SEXP Rfast2_fbed_reg(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_fedhc_skeleton(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_multinom_reg(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_weib_regs(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_welch_tests(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_wild_boot(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_mmhc_skeleton(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_negbin_reg(SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_negbin_regs(SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_gamma_regs(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP Rfast2_gamma_reg(SEXP, SEXP, SEXP, SEXP, SEXP); -//Stefanos +// Stefanos static const R_CallMethodDef CallEntries[] = { - {"Rfast2_benchmark", (DL_FUNC) &Rfast2_benchmark, 4}, - {"Rfast2_col_group", (DL_FUNC) &Rfast2_col_group, 4}, - {"Rfast2_colQuantile", (DL_FUNC) &Rfast2_colQuantile, 3}, - {"Rfast2_is_upper_tri", (DL_FUNC) &Rfast2_is_upper_tri, 2}, - {"Rfast2_is_lower_tri", (DL_FUNC) &Rfast2_is_lower_tri, 2}, - {"Rfast2_is_skew_symmetric", (DL_FUNC) &Rfast2_is_skew_symmetric, 1}, - {"Rfast2_lud", (DL_FUNC) &Rfast2_lud, 1}, - {"Rfast2_merge", (DL_FUNC) &Rfast2_merge, 2}, - {"Rfast2_Quantile", (DL_FUNC) &Rfast2_Quantile, 2}, - {"Rfast2_mmpc2", (DL_FUNC) &Rfast2_mmpc2, 10}, - {"Rfast2_rowQuantile", (DL_FUNC) &Rfast2_rowQuantile, 3}, - {"Rfast2_rowTrimMean", (DL_FUNC) &Rfast2_rowTrimMean, 3}, - {"Rfast2_colTrimMean", (DL_FUNC) &Rfast2_colTrimMean, 3}, - {"Rfast2_Runif", (DL_FUNC) &Rfast2_Runif, 3}, - {"Rfast2_Sample_int", (DL_FUNC) &Rfast2_Sample_int, 3}, - {"Rfast2_Sample", (DL_FUNC) &Rfast2_Sample, 3}, - {"Rfast2_trimmean", (DL_FUNC) &Rfast2_trimmean, 2}, + {"Rfast2_benchmark", (DL_FUNC)&Rfast2_benchmark, 4}, + {"Rfast2_col_group", (DL_FUNC)&Rfast2_col_group, 4}, + {"Rfast2_colQuantile", (DL_FUNC)&Rfast2_colQuantile, 4}, + {"Rfast2_is_upper_tri", (DL_FUNC)&Rfast2_is_upper_tri, 2}, + {"Rfast2_is_lower_tri", (DL_FUNC)&Rfast2_is_lower_tri, 2}, + {"Rfast2_is_skew_symmetric", (DL_FUNC)&Rfast2_is_skew_symmetric, 1}, + {"Rfast2_lud", (DL_FUNC)&Rfast2_lud, 1}, + {"Rfast2_merge", (DL_FUNC)&Rfast2_merge, 2}, + {"Rfast2_Quantile", (DL_FUNC)&Rfast2_Quantile, 3}, + {"Rfast2_mmpc2", (DL_FUNC)&Rfast2_mmpc2, 10}, + {"Rfast2_rowQuantile", (DL_FUNC)&Rfast2_rowQuantile, 4}, + {"Rfast2_rowTrimMean", (DL_FUNC)&Rfast2_rowTrimMean, 4}, + {"Rfast2_colTrimMean", (DL_FUNC)&Rfast2_colTrimMean, 4}, + {"Rfast2_Runif", (DL_FUNC)&Rfast2_Runif, 3}, + {"Rfast2_Sample_int", (DL_FUNC)&Rfast2_Sample_int, 3}, + {"Rfast2_Sample", (DL_FUNC)&Rfast2_Sample, 3}, + {"Rfast2_trimmean", (DL_FUNC)&Rfast2_trimmean, 3}, - {"Rfast2_inter", (DL_FUNC) &Rfast2_inter, 2}, - {"Rfast2_mmp_c", (DL_FUNC) &Rfast2_mmp_c, 10}, - - {"Rfast2_add_term", (DL_FUNC) &Rfast2_add_term, 9}, - {"Rfast2_colspml_mle", (DL_FUNC) &Rfast2_colspml_mle, 4}, - {"Rfast2_colcauchy_mle", (DL_FUNC) &Rfast2_colcauchy_mle, 4}, - {"Rfast2_colbeta_mle", (DL_FUNC) &Rfast2_colbeta_mle, 4}, - {"Rfast2_censweib_reg", (DL_FUNC) &Rfast2_censweib_reg, 5}, - {"Rfast2_fbed_reg", (DL_FUNC) &Rfast2_fbed_reg, 10}, - {"Rfast2_fedhc_skeleton", (DL_FUNC) &Rfast2_fedhc_skeleton, 7}, - {"Rfast2_multinom_reg", (DL_FUNC) &Rfast2_multinom_reg, 4}, - {"Rfast2_weib_regs", (DL_FUNC) &Rfast2_weib_regs, 6}, - {"Rfast2_welch_tests", (DL_FUNC) &Rfast2_welch_tests, 4}, - {"Rfast2_wild_boot", (DL_FUNC) &Rfast2_wild_boot, 7}, - {"Rfast2_mmhc_skeleton", (DL_FUNC) &Rfast2_mmhc_skeleton, 8}, - {"Rfast2_negbin_reg", (DL_FUNC) &Rfast2_negbin_reg, 4}, - {"Rfast2_negbin_regs", (DL_FUNC) &Rfast2_negbin_regs, 5}, - {"Rfast2_gamma_regs", (DL_FUNC) &Rfast2_gamma_regs, 6}, - {"Rfast2_gamma_reg", (DL_FUNC) &Rfast2_gamma_reg, 5}, - {NULL, NULL, 0} -}; + {"Rfast2_inter", (DL_FUNC)&Rfast2_inter, 2}, + {"Rfast2_mmp_c", (DL_FUNC)&Rfast2_mmp_c, 10}, + {"Rfast2_add_term", (DL_FUNC)&Rfast2_add_term, 9}, + {"Rfast2_colspml_mle", (DL_FUNC)&Rfast2_colspml_mle, 4}, + {"Rfast2_colcauchy_mle", (DL_FUNC)&Rfast2_colcauchy_mle, 4}, + {"Rfast2_colbeta_mle", (DL_FUNC)&Rfast2_colbeta_mle, 4}, + {"Rfast2_censweib_reg", (DL_FUNC)&Rfast2_censweib_reg, 5}, + {"Rfast2_fbed_reg", (DL_FUNC)&Rfast2_fbed_reg, 10}, + {"Rfast2_fedhc_skeleton", (DL_FUNC)&Rfast2_fedhc_skeleton, 7}, + {"Rfast2_multinom_reg", (DL_FUNC)&Rfast2_multinom_reg, 4}, + {"Rfast2_weib_regs", (DL_FUNC)&Rfast2_weib_regs, 6}, + {"Rfast2_welch_tests", (DL_FUNC)&Rfast2_welch_tests, 4}, + {"Rfast2_wild_boot", (DL_FUNC)&Rfast2_wild_boot, 7}, + {"Rfast2_mmhc_skeleton", (DL_FUNC)&Rfast2_mmhc_skeleton, 8}, + {"Rfast2_negbin_reg", (DL_FUNC)&Rfast2_negbin_reg, 4}, + {"Rfast2_negbin_regs", (DL_FUNC)&Rfast2_negbin_regs, 5}, + {"Rfast2_gamma_regs", (DL_FUNC)&Rfast2_gamma_regs, 6}, + {"Rfast2_gamma_reg", (DL_FUNC)&Rfast2_gamma_reg, 5}, + {NULL, NULL, 0}}; void R_init_Rfast2(DllInfo *info) { R_registerRoutines(info, NULL, CallEntries, NULL, NULL); R_useDynamicSymbols(info, FALSE); } - diff --git a/src/mmp_c.h b/src/mmp_c.h index 79fd601..73905b4 100644 --- a/src/mmp_c.h +++ b/src/mmp_c.h @@ -10,7 +10,7 @@ #include #include #include "cts.h" -#include "templates.h" +#include "Rfast2/templates.h" // [[Rcpp::plugins("cpp11")]] // [[Rcpp::depends("RcppArmadillo")]] diff --git a/src/mn.cpp b/src/mn.cpp index 07e0c22..500e58c 100644 --- a/src/mn.cpp +++ b/src/mn.cpp @@ -5,7 +5,7 @@ #include #include #include "mn.h" -#include "templates.h" +#include "Rfast2/templates.h" using namespace Rcpp; using namespace arma; diff --git a/src/mn.h b/src/mn.h index 91b4c46..2a50a67 100644 --- a/src/mn.h +++ b/src/mn.h @@ -2,7 +2,7 @@ // [[Rcpp::depends(RcppArmadillo)]] #include -#include "templates.h" +#include "Rfast2/templates.h" using namespace arma; using namespace Rcpp; diff --git a/src/multinom_reg.cpp b/src/multinom_reg.cpp index 5b7c28b..98e3b21 100644 --- a/src/multinom_reg.cpp +++ b/src/multinom_reg.cpp @@ -1,7 +1,7 @@ //Author: Stefanos Fafalios #include -#include "templates.h" +#include "Rfast2/templates.h" #include "apply_funcs_templates.h" #include #include "reg_lib_helper.h" diff --git a/src/negbin_reg.cpp b/src/negbin_reg.cpp index b471478..b4e67c3 100644 --- a/src/negbin_reg.cpp +++ b/src/negbin_reg.cpp @@ -5,7 +5,7 @@ #include "reg_lib2.h" #include "reg_lib_helper.h" #include "mn.h" -#include "templates.h" +#include "Rfast2/templates.h" using namespace Rcpp; using namespace arma; diff --git a/src/negbin_regs.cpp b/src/negbin_regs.cpp index c904bfb..be4dbde 100644 --- a/src/negbin_regs.cpp +++ b/src/negbin_regs.cpp @@ -5,7 +5,7 @@ #include "reg_lib2.h" #include "reg_lib_helper.h" #include "mn.h" -#include "templates.h" +#include "Rfast2/templates.h" using namespace Rcpp; using namespace arma; diff --git a/src/reg_lib2.cpp b/src/reg_lib2.cpp index a67a741..757ef13 100644 --- a/src/reg_lib2.cpp +++ b/src/reg_lib2.cpp @@ -3,7 +3,7 @@ #include #include "reg_lib2.h" #include -#include "templates.h" +#include "Rfast2/templates.h" #include "reg_lib_helper.h" #include "mn.h" #include "apply_funcs_templates.h" diff --git a/src/reg_lib_helper.cpp b/src/reg_lib_helper.cpp index ea576dc..e9ce725 100644 --- a/src/reg_lib_helper.cpp +++ b/src/reg_lib_helper.cpp @@ -3,7 +3,7 @@ #include #include "reg_lib_helper.h" #include -#include "templates.h" +#include "Rfast2/templates.h" #include "apply_funcs_templates.h" #include "reg_lib2.h" diff --git a/src/templates.h b/src/templates.h deleted file mode 100644 index 17d5b54..0000000 --- a/src/templates.h +++ /dev/null @@ -1,1188 +0,0 @@ -#include - -using namespace std; -using namespace arma; -using namespace Rcpp; - -#ifndef TEMPLATES_H -#define TEMPLATES_H - -//[[Rcpp::plugins(cpp11)]] - -template -struct pr{ - f first; - s second; - bool is_good; - pr(f first=0,s second=0):first(first),second(second),is_good(false){} -}; - -template -using Mfunction = T(*)(T,Args...); - -using Unary_Function=Mfunction; // typedef double (*Unary_Function)(double); // unary function -using Binary_Function=Mfunction; // typedef double (*Binary_Function)(double,double); // binary function -using Binary_Function_mat=Mfunction; // typedef double (*Binary_Function_mat)(mat,double); // binary function - - -/* - * F: unary function - * T1,T2: arguent class - * - * fill memory that starts from "startf" applying function F from "start" to "end" -*/ -template -void fill_with(T1 start,T1 end,T2 startf){ - for(;start!=end;++start,++startf){ - *startf=F(*start); - } -} - - -/* - * F: unary function - * T1,T2: arguent class - * - * fill memory that starts from "startf" applying function F to "x" -*/ -template -void fill_with(T1 x,T2 startf){ - for(typename T1::iterator start=x.begin();start!=x.end();++start){ - *startf=F(*start); - } -} - - -/* - * T: argument and return class - * - * applying function F to "x" and return "x". -*/ -template -T foreach(T x){ - for(typename T::iterator start=x.begin();start!=x.end();++start){ - *start=F(*start); - } - return x; -} - - -/* - * T1,T2: argument class - * - * count the frequency of "value" in "x" -*/ -template -int count_value_helper(T1 x,T2 value){ - int s=0; - for(typename T1::iterator start=x.begin();start!=x.end();++start){ - if(*start==value){ - ++s; - } - } - return s; -} - - -/* - * T: argument class - * - * count the frequency of "NAs" in "x" -*/ -template -int count_NAs(T x){ - int s=0; - for(typename T::iterator start=x.begin();start!=x.end();++start){ - if(R_IsNA(*start)){ - ++s; - } - } - return s; -} - - -/* - * T: argument class -*/ -template -double med_helper(typename T::iterator start,typename T::iterator last){ - double F; - const int sz=last-start,middle=sz/2-1; - if(sz%2==0){ - std::nth_element(start,start+middle,last); - F=(start[middle]+*(std::min_element(start+middle+1,last)))/2.0; - }else{ - std::nth_element(start,start+middle+1,last); - F=start[middle+1]; - } - return F; -} - - -/* - * Ret: return class - * T: argument class -*/ -template -Ret Order(T x,const bool stable,const bool descend,const int init_v){ - Ret ind(x.size()); - iota(ind.begin(),ind.end(),init_v); - if(descend){ - auto descend_func = [&](int i,int j){return x[i-init_v]>x[j-init_v];}; - stable ? stable_sort(ind.begin(),ind.end(),descend_func) : sort(ind.begin(),ind.end(),descend_func); - }else{ - auto func = [&](int i,int j){return x[i-init_v] -void min_max(T *start,T *end,T &mn, T &mx){ - T xxx; - mn=mx=*start; - start++; - for(;start!=end;++start){ - xxx=*start; - if(xxx>mx){ - mx=xxx; - } - else if(xxx -void maximum(T *start,T *end,T &mx){ - double xxx; - mx=*start; - start++; - for(;start!=end;++start){ - xxx=*start; - if(xxx>mx){ - mx=xxx; - } - } -} - - -/* - * T: argument class (not classes) -*/ -template -void minimum(T *start,T *end,T &mn){ - double xxx; - mn=*start; - start++; - for(;start!=end;++start){ - xxx=*start; - if(xxx -Ret Order_rank(T& x,const bool descend,const bool stable,const int n,const int k){ - Ret ind(x.size()-k); - iota(ind.begin(),ind.end(),0); - if(descend){ - auto descend_func = [&](int i,int j){return x[i]>x[j];}; - stable ? stable_sort(ind.begin(),ind.end()-n,descend_func) : sort(ind.begin(),ind.end()-n,descend_func); - }else{ - auto func = [&](int i,int j){return x[i] -double sum_with(T x){ - double a=0; - for(typename T::iterator start=x.begin();start!=x.end();++start){ - a+=F(*start); - } - return a; -} - -/* - * T: argument class - * - * applying function F to "x" and sum the elements of "x" -*/ -template F> -double sum_with(T *start,T *end){ - double a=0; - for(;start!=end;++start){ - a+=F(*start); - } - return a; -} - -/* - * F: a binary function - * T: argument class - * - * applying function F to "x" and sum the elements of "x" -*/ -template -double sum_with(T x,const double p){ - double a=0; - for(typename T::iterator start=x.begin();start!=x.end();++start){ - a+=F(*start,p); - } - return a; -} - -/* - * F: a binary function - * T: argument class - * - * applying function F to "x" and sum the elements of "x" -*/ -template -double sum_with(T1 x,T2& y){ - double a=0; - typename T1::iterator startx=x.begin(); - typename T2::iterator starty=y.begin(); - for(;startx!=x.end();++startx,++starty){ - a+=F(*startx,*starty); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T1 x,T2& y,Binary_Function F1,Binary_Function F2){ - double a=0; - typename T1::iterator startx=x.begin(); - typename T2::iterator starty=y.begin(); - for(;startx!=x.end();++startx,++starty){ - a=F2(a,F1(*startx,*starty)); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T1 x,T2& y){ - double a=0; - typename T1::iterator startx=x.begin(); - typename T2::iterator starty=y.begin(); - for(;startx!=x.end();++startx,++starty){ - a=F2(a,F1(*startx,*starty)); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T1 *x,T1 *endx,T2 *y){ - double a=0; - for(;x!=endx;++x,++y){ - a=F2(a,F1(*x,*y)); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T1 x,T2& y){ - double a=0; - typename T1::iterator startx=x.begin(); - typename T2::iterator starty=y.begin(); - for(;startx!=x.end();++startx,++starty){ - a=F3(a,F2(F1(*startx,*starty))); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T1 x,T2& y,Binary_Function F1,Unary_Function F2,Binary_Function F3){ - double a=0; - typename T1::iterator startx=x.begin(); - typename T2::iterator starty=y.begin(); - for(;startx!=x.end();++startx,++starty){ - a=F3(a,F2(F1(*startx,*starty))); - } - return a; -} - -/* - * F1: a binary function - * F2: a binary function - * T1: argument class - * T2: argument class - * - * applying function F1 to "x","y" and applying F2 to the result -*/ -template -double Apply(T x){ - double a=0; - typename T::iterator start=x.begin(); - for(;start!=x.end();++start){ - a=F2(a,F1(*start)); - } - return a; -} - -/* - * T: return and argument type (Rcpp and simple types. Not armadillo) -*/ -template -T square2(T x){ - return x*x; -} - -/* - * Ret: return class - * T: argument class -*/ -template -Ret Tabulate(T x,int &nroww){ - Ret f(nroww); - std::fill(f.begin(),f.end(),0); - typename Ret::iterator F=f.begin(); - typename T::iterator a=x.begin(); - for(;a!=x.end();++a){ - F[*a-1]++; - } - return f; -} - -/* - * Ret: return class - * T: argument class -*/ -template -Ret Tabulate(int* start,int* end,int& nroww){ - Ret f(nroww); - std::fill(f.begin(),f.end(),0); - typename Ret::iterator F=f.begin(); - for(;start!=end;++start){ - F[*start-1]++; - } - return f; -} - -/* - * Ret: return class - * T1: argument class - * T2: argument class -*/ -/* - * Ret: return class - * T1: argument class - * T2: argument class -*/ -template -Ret group_sum_helper(T1 x,T2 key,int *minn=nullptr,int *maxx=nullptr){ - int mn,mx; - const bool is_mn_null=(minn==nullptr),is_mx_null=(maxx==nullptr); - if(is_mx_null && is_mn_null){ - min_max(key.begin(),key.end(),mn,mx); - }else if(is_mx_null){ - mn=*minn; - maximum(key.begin(),key.end(),mx); - }else if(is_mn_null){ - mx=*maxx; - minimum(key.begin(),key.end(),mn); - }else{ - mn=*minn; - mx=*maxx; - } - typename T2::iterator kk=key.begin(); - vector f(mx-mn+1); - vector is_good(mx-mn+1); - typename T1::iterator xx=x.begin(),rr; - vector::iterator ff=f.begin(); - vector::iterator ok; - int index; - for(;xx!=x.end();++xx,++kk){ - index = *kk-mn; - is_good[index]=true; - f[index] += *xx; - } - int number_of_values=0; - for(auto v : is_good){ - if(v) - ++number_of_values; - } - Ret res(number_of_values); - for(rr=res.begin(),ok=is_good.begin(),ff=f.begin();ff!=f.end();++ff){ - if(*ok++){ - *rr++=*ff; - } - } - return res; -} - - -/* Distixos gia tora to T prepei na einai panta mat - * dioti mono ta class tou arma exoun n_rows,n_cols - * Ret: return class - * T: argument class -*/ -template -Ret cross_x_y(T x,T y){ - const int ncl=x.n_cols,nrw=x.n_rows,p=y.n_cols; - Ret f(ncl,p); - Vec yv(nrw); - int i,j; - for(i=0;i -Ret cross_x(T x){ - const int ncl=x.n_cols; - Ret f(ncl,ncl); - double a; - int i,j; - for(i=0;i -Ret design_matrix_helper(T x) { - int i=0; - const int n=x.size(); - T tmp=sort_unique(x); - typename T::iterator xx=x.begin(),leksi_bg,leksi_en; - Ret Final(n,tmp.size()); - std::fill(Final.begin(),Final.end(),0); - for(leksi_bg=tmp.begin(),leksi_en=tmp.end(),i=0;xx!=x.end();++xx,++i){ - Final(i,lower_bound(leksi_bg,leksi_en,*xx)-leksi_bg)=1; - } - return Final; -} - -template -inline T madd(T x,T y){ - return x+y; -} - - -template -inline T mless(T x,T y){ - return x -inline T mless_eq(T x,T y){ - return x<=y; -} - -template -inline T mgreater_eq(T x,T y){ - return x>=y; -} - -template -inline T mgreater(T x,T y){ - return x>y; -} - - -template -inline T mdiv(T x,T y){ - return x/y; -} - - -template -inline T mdiff(T x,T y){ - return x-y; -} - - -template -inline T mmult(T x,T y){ - return x*y; -} - -template -inline T mequal(T x,T y){ - return x==y; -} - -template -inline T mmax(T x,T y){ - return std::max(x,y); -} - -template -inline T mmin(T x,T y){ - return std::min(x,y); -} - -template -inline T mand(T x,T y){ - return x&&y; -} - -template -inline T mor(T x,T y){ - return x||y; -} - -template -inline RET madd(T1 x,T2 y){ - return x+y; -} - - -template -inline RET mless(T1 x,T2 y){ - return x -inline RET mless_eq(T1 x,T2 y){ - return x<=y; -} - -template -inline RET mgreater_eq(T1 x,T2 y){ - return x>=y; -} - -template -inline RET mgreater(T1 x,T2 y){ - return x>y; -} - - -template -inline RET mdiv(T1 x,T2 y){ - return x/y; -} - - -template -inline RET mdiff(T1 x,T2 y){ - return x-y; -} - - -template -inline RET mmult(T1 x,T2 y){ - return x*y; -} - -template -inline RET mequal(T1 x,T2 y){ - return x==y; -} - -template -inline RET mmax(T1 x,T2 y){ - return std::max(x,y); -} - -template -inline RET mmin(T1 x,T2 y){ - return std::min(x,y); -} - -template -inline RET mand(T1 x,T2 y){ - return x&&y; -} - -template -inline RET mor(T1 x,T2 y){ - return x||y; -} - - -template -void myoperator(T f[],T &x,T *y,int &len){ - T *ff=f; - for(int i=0;i -void as_integer_h_sorted(vector x,IntegerVector &f,const int init,const T val){ - const int n=x.size(); - int i,j=0,c=init; - sort(x.begin(),x.end()); - auto v=x[j]; - f[0]=init; - for(i=1;i -void as_integer_h(vector x,IntegerVector &f,const int init,const T val){ - const int n=x.size(); - int i,j=0,c=init; - vector ind=Order< vector,vector >(x,false,false,0); // diorthoseiii - x.push_back(val); - T v=x[ind[j]]; - f[ind[0]]=init; - for(i=1;i -void as_integer_h_with_names(vector x,List &L,const int init,const T val){ - const int n=x.size()+1; - int i,j=0,c=init; - vector ind=Order< vector,vector >(x,false,false,0); // diorthoseiii - x.push_back(val); - ind.push_back(0); - vector w; - T v=x[ind[j]]; - IntegerVector f(n-1); - f[ind[0]]=init; - for(i=1;i -vector table_use_na(vector x,const int use_na){ - typename vector::iterator new_end = std::remove_if(x.begin(),x.end(),R_IsNA); - sort(x.begin(),new_end); - vector n; - typename vector::iterator a=x.begin(),b=a+1; - List l; - int s=1; - for(;b!=new_end;++b){ - if(*a!=*b){ - n.push_back(s); - a=b; - s=1; - }else - ++s; - } - if(use_na==1){ - n.push_back(new_end-x.begin()); - } - return n; -} - -template -vector table_simple(vector x){ - sort(x.begin(),x.end()); - x.push_back(T(0)); - vector n; - typename vector::iterator a=x.begin(),b=a+1; - int s=1; - for(;b!=x.end();++b){ - if(*a!=*b){ - n.push_back(s); - a=b; - s=1; - }else - ++s; - } - return n; -} - -template -void table2_like_r(vector x,vector y,IntegerMatrix &f,const T val){ - //cout<<"table2_like_r\n"; - const int n=x.size(); - int mx_x,mx_y; - IntegerVector ix(n),iy(n); - as_integer_h(x,ix,0,val); - as_integer_h(y,iy,0,val); - maximum(ix.begin(),ix.end(),mx_x); - maximum(iy.begin(),iy.end(),mx_y); - f=IntegerMatrix(mx_x+1,mx_y+1); - for(int i=0;i -void table2_like_r_with_names(vector x,vector y,List &L,const T val){ - //cout<<"table2_like_r_with_names"<(x,lx,0,val); - as_integer_h_with_names(y,ly,0,val); - IntegerVector ix=lx["f"],iy=ly["f"]; - maximum(ix.begin(),ix.end(),mx_x); - maximum(iy.begin(),iy.end(),mx_y); - IntegerMatrix f(mx_x+1,mx_y+1); - for(int i=0;i -Ret rank_mean(T x,const bool descend){ - const int n=x.size(),n_1=n+1; - int i,j=0; - x.resize(n_1); - x[n]=0; - I ind=Order_rank(x,descend,false,1,0); - int k=0,m,times=0; - Ret f(n); - double mn=0.0,v=x[ind[j]]; - for(i=1;i -Ret rank_max(T x,const bool descend){ - const int n=x.size(),n_1=n+1; - int i,j=0; - x.resize(n_1); - x[n]=0; - I ind=Order_rank(x,descend,false,1,0); - int k=0,m,times=0; - Ret f(n); - double v=x[ind[j]]; - for(i=1;i -Ret rank_min(T x,const bool descend){ - const int n=x.size(); - int i,j=0; - I ind=Order_rank(x,descend,false,0,1); - Ret f(n); - double v=x[ind[j]]; - f[ind[0]]=1; - for(i=1;i -Ret rank_first(T x,const bool descend,const bool stable){ - const int n=x.n_elem; - I ind=Order_rank(x,descend,stable,0,1); - Ret f(n); - for(int i=0;i -NumericVector eachcol_apply_helper(NumericMatrix& x,NumericVector& y,SEXP ind){ - const bool is_ind_null = Rf_isNull(ind); - const int n = is_ind_null ? x.ncol() : LENGTH(ind); - NumericVector f(n); - if(is_ind_null){ - for(int i=0;i::Column,NumericVector,oper,func >(x.column(i),y); - } - }else{ - IntegerVector indd(ind); - for(int i=0;i::Column,NumericVector,oper,func >(x.column(indd[i]-1),y); - } - } - return f; -} - -template -double apply_eachrow_helper(SEXP x,SEXP y){ - int ncol=Rf_ncols(x),nrow=Rf_nrows(x); - SEXP mat=Rf_duplicate(x); - double *xx=REAL(mat),*end=xx+ncol*nrow,*yy=REAL(y),y3,*x3,s=0; - for(;xx!=end;++yy){ - y3=*yy; - for(x3=xx,xx+=nrow;x3!=xx;++x3){ - s=func(s,oper(*x3,y3)); - } - } - return s; -} - -template -SEXP eachrow_helper(SEXP x,SEXP y){ - int ncol=Rf_ncols(x),nrow=Rf_nrows(x); - SEXP mat=PROTECT(Rf_allocMatrix(type,nrow,ncol)); - T *xx=(T *) DATAPTR(x),*xend=xx+ncol*nrow,*yy=(T *) DATAPTR(y),yvalue,*x3; - RETURN_TYPE *m=(RETURN_TYPE*)DATAPTR(mat); - for(;xx!=xend;++yy){ - yvalue=*yy; - for(x3=xx,xx+=nrow;x3!=xx;++x3,++m){ - *m=oper(*x3,yvalue); - } - } - UNPROTECT(1); - return mat; -} - -template -double sum_x_op_y(SEXP x,SEXP y){ - const int n=LENGTH(x); - double s=0,*startx=REAL(x),*endx=startx+n,*starty=REAL(y); - for(;startx!=endx;++startx,++starty){ - s=func(s,oper(*startx,*starty)); - } - return s; -} - -template -double sum_x_op_x(SEXP x){ - const int n=LENGTH(x); - double s=0,*startx=REAL(x),*endx=startx+n; - for(;startx!=endx;++startx){ - s=func(s,oper(*startx,*startx)); - } - return s; -} - -//sort_unique,len_sort_unique,sort_int -template -void max_neg_pos(T *start,T *end,T &mx,T &mn,bool &pos,bool &neg){ - mn=mx=*start; - T x; - for(;start!=end;++start){ - x=*start; - if(x<0){ - neg=true; - if(mn>x){ - mn=x; - } - }else{ - pos=true; - if(mx Cond,Mfunction Oper> -int Apply_helper(int *start,int *end,int& val){ - int t=0; - for(;start!=end;++start){ - if(Cond(*start,val)){ - t=Oper(t,*start); - } - } - return t; -} - -template // less or greater -NumericVector negative_or_positive(NumericVector &x){ - double v,val=x[0]; - for(auto xx=x.begin()+1;xx!=x.end();++xx){ - v=*xx; - if(less_or_greater(v,0)){ - if(min_or_max(v,val)){ - val=v; - } - } - } - return NumericVector::create(val); -} - -template // less or greater -NumericVector negative_or_positive_min_max(NumericVector &x){ - double v,mn=x[0],mx=mn; - for(auto xx=x.begin()+1;xx!=x.end();++xx){ - v=*xx; - if(less_or_greater(v,0)){ - if(vmx){ - mx=v; - } - } - } - return NumericVector::create(mn,mx); -} - -template -double nth_simple(T& x,const int elem,const bool descend){ - descend ? - nth_element(x.begin(),x.begin()+elem-1,x.end(),[&](double a,double b){return a>b;}) - : - nth_element(x.begin(),x.begin()+elem-1,x.end()); - - return x[elem-1]; -} - -template -double nth_na_rm(T& x,const int& elem,const bool& descend){ - const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); - descend ? - nth_element(x.begin(),x.begin()+((elemb;}) - : - nth_element(x.begin(),x.begin()+((elem -double nth_helper(T& x,const int elem,const bool descend,const bool na_rm){ - return na_rm ? nth_na_rm(x,elem,descend) : nth_simple(x,elem,descend); -} - - -template -int nth_index_simple(T& x,const int& elem,const bool& descend){ - IntegerVector ind=seq(1,x.size()); - descend ? - nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1]>x[j-1];}) - : - nth_element(ind.begin(),ind.begin()+elem-1,ind.end(),[&](int i,int j){return x[i-1] -int nth_index_na_rm(T& x,const int& elem,const bool& descend){ - const int new_end=remove_if(x.begin(),x.end(),R_IsNA)-x.begin(); - IntegerVector ind= seq(1,new_end); - descend ? - nth_element(ind.begin(),ind.begin()+((elemx[j-1];}) - : - nth_element(ind.begin(),ind.begin()+((elem -int nth_helper_index(T& x,const int elem,const bool descend,const bool na_rm){ - return na_rm ? nth_index_na_rm(x,elem,descend) : nth_index_simple(x,elem,descend); -} - -template func,const int init_val=0> -SEXP group_col_h(SEXP x,SEXP gr,const int length_unique){ - const int ncl=Rf_ncols(x),nrw=Rf_nrows(x); - SEXP f=PROTECT(Rf_allocMatrix(TYPEOF(x),length_unique,ncl)); - int *ggr=INTEGER(gr); - T *ff=(T*)DATAPTR(f),*xx=(T*)DATAPTR(x); - for(int j=0;j -SEXP group_col_med_h(SEXP x,SEXP gr,const int length_unique){ - const int ncl=Rf_ncols(x),nrw=Rf_nrows(x); - SEXP f=PROTECT(Rf_allocMatrix(TYPEOF(x),length_unique,ncl)); - int *ggr=INTEGER(gr); - T *ff=(T*)DATAPTR(f),*xx=(T*)DATAPTR(x); - vector> eachcol_mat(length_unique,vector()); - for(int j=0;j& tmp = eachcol_mat[i]; - ff[i+col_index_f]=med_helper>(tmp.begin(),tmp.end()); - tmp.clear(); - } - } - UNPROTECT(1); - return f; -} - -template -Ret Quantile(T x,colvec& probs){ - Ret f(probs.n_elem); - - if(probs.size() > 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 0.5){ - a = nth_simple(x, hf-1,false); - b = *min_element(x.begin()+hf,x.end()); - }else{ - b = nth_simple(x, hf,false); - a = *max_element(x.begin(),x.begin()+hf); - } - f[i] = a + (h-hf) * ( b - a ); - } - } - - return f; -} - -template -double trimmean_h(T x,const double a) { - const int n = x.n_elem; - int b1 = a * n; - int b11 = std::ceil(b1); - if ( b1 - b11 == 0 ) { - b1 = b11 + 1; - } else{ - b1 = b11; - } - const double a1 = nth_simple(x, b1,false); - const double a2 = nth_simple(x, n - b1 + 1,false); - double s=0; - int p=0; - for(auto xx : x){ - if(xx>=a1 && xx<=a2){ - s+=xx; - ++p; - } - } - return s/p; -} - - -#endif - - diff --git a/src/utilities.cpp b/src/utilities.cpp index 90dfb72..d3007c2 100644 --- a/src/utilities.cpp +++ b/src/utilities.cpp @@ -4,20 +4,15 @@ // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 #include -#include #include #include -#include "templates.h" +#include "Rfast2.h" using namespace Rcpp; using namespace arma; -using std::sort; -using std::vector; using std::merge; using std::vector; -using std::min; - List lud(NumericMatrix x){ int ncl=x.ncol(),nrw=x.nrow(),i=0,j; @@ -149,35 +144,30 @@ END_RCPP ////////////////////////////////////////////////////////////////////////////////////////////////////////// -NumericVector Quantile(colvec x,NumericVector Probs){ - colvec probs(Probs.begin(),Probs.size(),false); - return Quantile(x,probs); -} - -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< colvec >::type x(xSEXP); - traits::input_parameter< NumericVector >::type Probs(ProbsSEXP); - __result = Quantile(x,Probs); + traits::input_parameter< const bool >::type parallel(parallelSEXP); + NumericVector x(xSEXP); + NumericVector Probs(ProbsSEXP); + + colvec probs(Probs.begin(),Probs.size(),false); + __result = Rfast::Quantile(x,probs,parallel); return __result; END_RCPP } /***********************************************************************************/ -double trimmean(colvec x,const double a){ - return trimmean_h(x,a); -} - -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 = trimmean(x,a); + traits::input_parameter< const bool >::type parallel(parallelSEXP); + __result = Rfast::TrimMean(x,a,parallel); return __result; END_RCPP } diff --git a/src/weib_regs.cpp b/src/weib_regs.cpp index 38c0a2c..61fd526 100644 --- a/src/weib_regs.cpp +++ b/src/weib_regs.cpp @@ -4,7 +4,7 @@ #include #include "reg_lib2.h" #include "reg_lib_helper.h" -#include "templates.h" +#include "Rfast2/templates.h" #include using namespace Rcpp; diff --git a/src/welch_tests.cpp b/src/welch_tests.cpp index 8c98cd0..705fb58 100644 --- a/src/welch_tests.cpp +++ b/src/welch_tests.cpp @@ -2,7 +2,7 @@ #include #include -#include "templates.h" +#include "Rfast2/templates.h" // [[Rcpp::depends(RcppArmadillo)]] using namespace arma; diff --git a/src/wild_boot.cpp b/src/wild_boot.cpp index 6d76219..e927eb5 100644 --- a/src/wild_boot.cpp +++ b/src/wild_boot.cpp @@ -2,7 +2,7 @@ #include #include -#include "templates.h" +#include "Rfast2/templates.h" #include "mn.h" using namespace Rcpp;