From 8d17029e779f0c2f99f5161f23114026c319c8d3 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 21:48:46 +0100 Subject: [PATCH 01/14] Add constant mpolys handling in -.mpoly This may be unnecessary as it's done in +.mpoly --- R/mpolyArithmetic.R | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/R/mpolyArithmetic.R b/R/mpolyArithmetic.R index 9b4ce29..dacc6c0 100644 --- a/R/mpolyArithmetic.R +++ b/R/mpolyArithmetic.R @@ -63,14 +63,25 @@ NULL #' @rdname mpolyArithmetic #' @export `-.mpoly` <- function(e1, e2){ - + + ## if either is constant, mpoly it + if(!is.list(e1)){ + stopifnot(is.numeric(e1) && length(e1) == 1) + e1 <- list(c(coef = e1)) + } + + if(!is.list(e2)){ + stopifnot(is.numeric(e2) && length(e2) == 1) + e2 <- list(c(coef = e2)) + } + ## flip coefficients of each term e2 <- lapply(e2, function(v){ v["coef"] <- -v["coef"] v }) class(e2) <- "mpoly" - + ## add e1 + e2 } From 85e027b2070ec56ffe5c9d092f0d28e3dceb1342 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 21:50:37 +0100 Subject: [PATCH 02/14] Add Rcpp dependency to DESCRIPTION --- DESCRIPTION | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ecf56b7..5cc95fc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,8 +4,10 @@ Title: Symbolic Computation and More with Multivariate Polynomials Version: 1.0.5.902 URL: https://github.com/dkahle/mpoly BugReports: https://github.com/dkahle/mpoly/issues -Authors@R: person("David", "Kahle", email = "david.kahle@gmail.com", role = c("aut", "cre")) +Authors@R: c(person("David", "Kahle", email = "david.kahle@gmail.com", role = c("aut", "cre")), + person("Blagoje", "Ivanovic", email = "blagoje.r.ivanovic@gmail.com", role = c("aut"))) Description: Symbolic computing with multivariate polynomials in R. +LinkingTo: Rcpp Imports: stringr (>= 1.0.0), partitions, @@ -14,10 +16,12 @@ Imports: ggplot2, polynom, orthopolynom, - tidyr + tidyr, + Rcpp Suggests: testthat, magrittr, dplyr +SystemRequirements: C++11 License: GPL-2 RoxygenNote: 6.0.1 From 3f8e8f3991430b02cb2a31d014d18fe6a59103be Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 21:50:49 +0100 Subject: [PATCH 03/14] Add Rcpp dependency to help.R --- R/help.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/help.R b/R/help.R index 23a9053..0708e75 100644 --- a/R/help.R +++ b/R/help.R @@ -1,7 +1,9 @@ #' Multivariate polynomials in R. #' #' A package for symbolic computation and more with multivariate polynomials -#' +#' @import Rcpp +#' @importFrom Rcpp sourceCpp +#' @useDynLib mpoly #' @importFrom stringr str_c str_detect str_sub str_replace_all str_replace #' str_count str_dup str_extract str_extract_all str_locate str_locate_all #' str_split str_sub<- str_trim fixed From 252efe80cee8c26c2e79caec03a609842ef7a6e7 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 21:58:05 +0100 Subject: [PATCH 04/14] Add helper functions for cpp code --- inst/include/helpers.h | 24 ++++++++++++++++++++++++ src/helpers.cpp | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+) create mode 100644 inst/include/helpers.h create mode 100644 src/helpers.cpp diff --git a/inst/include/helpers.h b/inst/include/helpers.h new file mode 100644 index 0000000..93b4e68 --- /dev/null +++ b/inst/include/helpers.h @@ -0,0 +1,24 @@ +#ifndef mpoly_ADD_H +#define mpoly_ADD_H + +using namespace Rcpp; + +// concatenate string vectors +std::vector scat(std::vector x, + std::vector y); + +// concatenate two numeric vectors (as c() in R) +NumericVector vcat(const NumericVector& x, const NumericVector& y); + +// helper function to convert types to string +template +std::string to_string(T const& value) { + std::stringstream sstr; + sstr << value; + return sstr.str(); +}; + +// Vectorized equality operator for string vs Character vector +LogicalVector vecstreq(const CharacterVector& vec, std::string str); + +#endif diff --git a/src/helpers.cpp b/src/helpers.cpp new file mode 100644 index 0000000..b279c38 --- /dev/null +++ b/src/helpers.cpp @@ -0,0 +1,35 @@ +#include +#include +using namespace Rcpp; + +// concatenate string vectors +std::vector scat(std::vector x, + std::vector y) { + std::vector z(x.size() + y.size()); + std::copy(x.begin(), x.end(), z.begin()); + std::copy(y.begin(), y.end(), z.begin() + x.size()); + return z; +} + +// concatenate two numeric vectors (as c() in R) +NumericVector vcat(const NumericVector& x, const NumericVector& y) { + NumericVector z(x.size() + y.size()); + + std::copy(x.begin(), x.end(), z.begin()); + std::copy(y.begin(), y.end(), z.begin() + x.size()); + + z.names() = scat(x.names(), y.names()); + + return z; +} + +// Vectorized equality operator for string vs Character vector +LogicalVector vecstreq(const CharacterVector& vec, std::string str) { + LogicalVector res(vec.size()); + + for (int i = 0; i < vec.size(); i++) { + res[i] = vec[i] == str; + } + + return res; +} From e94e49b69894f59f4f729f3a8b2354975196f0d8 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 21:58:13 +0100 Subject: [PATCH 05/14] Add Makevars to include inst/include folder during compilation --- src/Makevars | 1 + 1 file changed, 1 insertion(+) create mode 100644 src/Makevars diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..91dc105 --- /dev/null +++ b/src/Makevars @@ -0,0 +1 @@ +PKG_CPPFLAGS = -I../inst/include/ \ No newline at end of file From a45064536b5e832e2d5098f5342edfb686b0e0b1 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:00:01 +0100 Subject: [PATCH 06/14] Add mpoly.cpp This file contains the cpp function reducePoly which handles the bulk of the mpoly function in R, speeding it up quite a bit. --- src/mpoly.cpp | 217 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 217 insertions(+) create mode 100644 src/mpoly.cpp diff --git a/src/mpoly.cpp b/src/mpoly.cpp new file mode 100644 index 0000000..084d2d3 --- /dev/null +++ b/src/mpoly.cpp @@ -0,0 +1,217 @@ +#include +#include +#include +#include +using namespace Rcpp; + +// [[Rcpp::plugins(cpp11)]] + + +NumericVector combineDegrees(const NumericVector& vec) { + CharacterVector names = vec.names(); + CharacterVector unames = unique(names); + NumericVector combined(unames.size(), 0); + + combined.names() = unames; + combined["coef"] = 1; + + for (int i = 0; i < vec.size(); i++) { + if (names[i] == "coef") { + combined["coef"] = combined["coef"] * vec[i]; + } + else { + combined[(std::string)names[i]] = combined[(std::string)names[i]] + vec[i]; + } + } + + // remove zero degree elements + NumericVector res = combined[ combined != 0 | vecstreq(unames, "coef") ]; + + return res; +} + + +// Variable comparison function. A variable is before another if it comes first in vars vector. +// coef is after (not before) everything +struct compareVarOrder { + compareVarOrder(const CharacterVector& vars) { this->vars = vars; } + + bool operator()(std::string var1, std::string var2) const + { + if (var1 == "coef") return false; + if (var2 == "coef") return true; + + for (int i = 0; i < vars.size(); i++) { + if (vars[i] == var1) return true; + if (vars[i] == var2) return false; + } + return true; + } + + CharacterVector vars; +}; + + +// [[Rcpp::export]] +NumericVector sortVars(const NumericVector& term, const CharacterVector& vars) { + std::vector names = term.names(); + std::sort(names.begin(), names.end(), compareVarOrder(vars)); + + NumericVector sorted(term.size()); + for (int i = 0; i < term.size(); i++) { + sorted[i] = term[names[i]]; + } + + sorted.names() = names; + + return sorted; +} + + +// [[Rcpp::export]] +CharacterVector getMonomialStrings(const List& poly) { + CharacterVector res(poly.size()); + + for (int i = 0; i < poly.size(); i++) { + NumericVector term = poly[i]; + CharacterVector names = term.names(); + std::string monom = ""; + + if (term.size() == 1 && names[0] == "coef") { + monom = "coef"; + } + else { + for (int j = 0; j < term.size(); j++) { + if (!(names[j] == "coef")) + monom += (std::string)names[j] + to_string(term[j]); + } + } + + res[i] = monom; + } + + return res; +} + + +// [[Rcpp::export]] +List combineMonomials(const List& poly) { + CharacterVector monomials = getMonomialStrings(poly); + CharacterVector uniqueMonomials = unique(monomials); + int nUnique = uniqueMonomials.size(); + int nMonom = monomials.size(); + + if (nMonom == nUnique) + return poly; + + // make a map of (unique) monomial strings and their coefficients + std::unordered_map monomCoefMap; + // and a monomial string to index in monomials vector + // we need tese indices to be able to pull vector from poly to use as a template + // and only change its coefficient + std::unordered_map monomIndexMap; + + // populate the maps with summed coefficient values and indices + for (int i = 0; i < nMonom; i++) { + NumericVector term = poly[i]; + monomCoefMap[(std::string)monomials[i]] += term["coef"]; + monomIndexMap[(std::string)monomials[i]] = i; + } + + // resulting list + List res(nUnique); + + // populate resulting list with summed coefficients + for (int i = 0; i < nUnique; i++) { + NumericVector templateTerm = poly[ monomIndexMap[(std::string)uniqueMonomials[i]] ]; + templateTerm["coef"] = monomCoefMap[(std::string)uniqueMonomials[i]]; + res[i] = templateTerm; + } + + return res; +} + +bool isSubset(CharacterVector sub, CharacterVector sup) { + for (int i = 0; i < sub.size(); i++) { + if(std::find(sup.begin(), sup.end(), sub[i]) == sup.end()) + return false; + } + return true; +} + +// [[Rcpp::export]] +List tidyTerms(const List& poly, const CharacterVector& vars) { + List res(poly.size()); + + // static vector which we'll use to add coef to terms that don't have it + NumericVector coef1vec = NumericVector::create(Named("coef", 1.0)); + + for (int i = 0; i < poly.size(); i++) { + NumericVector term = poly[i]; + CharacterVector names = term.names(); + + // add coef = 1 if no coef present, but vars are + if (!term.containsElementNamed("coef") && isSubset(names, vars)) + term = vcat(term, coef1vec); + + // merge same variables by adding up their degrees + term = combineDegrees(term); + + // sort the variables + term = sortVars(term, vars); + + res[i] = term; + } + + return res; +} + +// [[Rcpp::export]] +List removeZeros(const List& poly){ + // to avoid expanding the list, we'll first check how many zeros there are + // and then proceed to populate the resulting list without zero elements + int nZeros = 0; + + for (int i = 0; i < poly.size(); i++) { + NumericVector term = poly[i]; + if (term["coef"] == 0.0) { + nZeros++; + } + } + + if (nZeros == 0) + return poly; + + List res(poly.size() - nZeros); + + int j = 0; + for (int i = 0; i < poly.size(); i++) { + NumericVector term = poly[i]; + if (term["coef"] != 0) { + res[j] = term; + j++; + } + } + + return res; +} + + +// [[Rcpp::export]] +List reducePoly(const List& poly, const CharacterVector& vars){ + // tidy up the terms in the poly list + List tidyPoly = tidyTerms(poly, vars); + + // combine the monomials in the polynomial, i.e. add up same monomials + List combinedPoly = combineMonomials(tidyPoly); + + // remove elements with zero coefficients + List reducedPoly = removeZeros(combinedPoly); + + // if everything is removed, we have a zero polynomial and thus return that + if (reducedPoly.size() == 0) { + reducedPoly = List::create(NumericVector::create(Named("coef", 0.0))); + } + + return reducedPoly; +} From 7795ec52f1dfebd69685ee894f2852dc818d3a37 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:00:27 +0100 Subject: [PATCH 07/14] Make mpoly.R use the Rcpp implementation --- R/mpoly.R | 153 +----------------------------------------------------- 1 file changed, 1 insertion(+), 152 deletions(-) diff --git a/R/mpoly.R b/R/mpoly.R index 21a4d54..7e47e6d 100644 --- a/R/mpoly.R +++ b/R/mpoly.R @@ -50,42 +50,7 @@ mpoly <- function(list, varorder){ flatList <- flatList[names(flatList) != "coef"] if(!all(is.wholenumber(flatList)) || any(flatList < 0)){ stop("degrees must be nonnegative integers.", call. = FALSE) - } - - - # give terms without a coefficient a coef of 1 - addCoefIfMissing <- function(v){ - if(any("coef" == names(v))) return(v) - c(v, coef = 1) } - list <- lapply(list, addCoefIfMissing) - - - - ## organize - - # remove terms with coef 0 - list <- filterOutZeroTerms(list) - - - # remove 0 degrees, combine like degrees, single coef as rightmost element - list <- lapply(list, function(v){ - # separate vardegs from coefs - coef_ndx <- which(names(v) == "coef") - coefs <- v[coef_ndx] - v <- v[-coef_ndx] - - # combine like degrees (sum) - if(length(names(v)) != length(unique(names(v)))) v <- fastNamedVecTapply(v, sum) - - # combine like coefficients (product) - coefs <- c(coef = prod(coefs)) - - # remove zero degree elements, combine and return - c(v[v != 0], coefs) - }) - - ## set intrinsic varorder - done again after 0 degrees are removed vars <- unique(names(flatList)) @@ -101,126 +66,10 @@ mpoly <- function(list, varorder){ } vars <- varorder } - - - # sort variables in terms - list <- lapply(list, function(v){ - p <- length(v) - 1L - if(p == 0L) return(v) - c( (v[1:p])[intersect(vars, names(v[1:p]))], v["coef"] ) - }) - - - - ## prepare to check if like terms are present - monomials <- vapply(list, function(v){ - p <- length(v) - 1 # remove coef on monomials - paste(names(v[1:p]), v[1:p], sep = "", collapse = "") - }, character(1)) - # e.g. c("x1", "y1", "y1", "x2z1", "y4t3", "x1", "x1", "coef5") - - unique_monomials <- unique(monomials) - - - - ## check if like terms are present and, if so, correct - if(length(monomials) != length(unique_monomials)){ - - matchedMonomials <- match(monomials, unique_monomials) - matchedMonomials <- factor(matchedMonomials, levels = 1:max(matchedMonomials)) - ndcs2combine <- split.default(1:length(list), matchedMonomials) - - list <- lapply(ndcs2combine, function(v){ - if(length(v) == 1) return(list[[v]]) - coef <- sum(vapply(list[v], `[`, double(1), length(list[[v[1]]]))) - v <- list[[v[1]]] - v["coef"] <- coef - v - }) - - names(list) <- NULL # i.e. list <- unname(list) - } - - - - ## combine constant terms - ## mpoly(list(c(x = 1, coef = 1), c(coef = 1), c(coef = 2))) - nonConstantTerms <- fastFilter(isNotLengthOne, list) - constantTerms <- fastFilter(isLengthOne, list) - - if(length(constantTerms) > 0){ - list <- c( - nonConstantTerms, - list(Reduce(`+`, constantTerms)) - ) - } - - - - ## re-organize after like-terms combined - list <- filterOutZeroTerms(list) - + list <- reducePoly(list, vars) ## return classed list class(list) <- "mpoly" list } - - - - - - - - - - - -filterOutZeroTerms <- function(list){ - - # a term is zero if any of its coefficients are 0 - hasNonZeroCoef <- function(v) all(v["coef"] != 0) - - # remove terms with coef 0 - list <- fastFilter(hasNonZeroCoef, list) - - # if all terms have been eliminated, recreate it with a 0 coef - if(length(list) == 0) list <- list(c(coef = 0)) - - # return - list -} - - - - - - - -isLengthOne <- function(x) length(x) == 1L -isNotLengthOne <- function(x) length(x) > 1L -fastFilter <- function(f, x) x[vapply(x, f, logical(1))] - - - - - - - -fastNamedVecTapply <- function(x, f, type = double(1)){ - uniqueNames <- unique(names(x)) - matchedNames <- match(names(x), uniqueNames) # indices - matchedNames <- factor(matchedNames, levels = 1:max(matchedNames)) - groupIndices <- split.default(1:length(x), matchedNames) - out <- vapply(groupIndices, function(ndcs) f(x[ndcs]), type) - names(out) <- uniqueNames - out -} -# x <- 1:10 -# names(x) <- sample(letters[1:3], 10, replace = TRUE) -# fastNamedVecTapply(x, sum) -# tapply(x, names(x), sum) - - - From 4cb21208efa7b035c1d65601426918eecb008e87 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:01:24 +0100 Subject: [PATCH 08/14] Add all.equal.mpoly function This makes tests pass for mpoly, as it makes all.equal not compare lists, but use the equality of polynomials --- R/mpolyEqual.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/mpolyEqual.R b/R/mpolyEqual.R index 8df0fd2..2089d2b 100644 --- a/R/mpolyEqual.R +++ b/R/mpolyEqual.R @@ -51,5 +51,7 @@ } - +all.equal.mpoly <- function(e1, e2) { + e1 == e2 +} From 4fa0b71eca1049aa4318178da40bb1637344fa92 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:02:22 +0100 Subject: [PATCH 09/14] Add Rcpp versions of mpoly arithmetic --- src/mpolyArithmetic.cpp | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 src/mpolyArithmetic.cpp diff --git a/src/mpolyArithmetic.cpp b/src/mpolyArithmetic.cpp new file mode 100644 index 0000000..21cd445 --- /dev/null +++ b/src/mpolyArithmetic.cpp @@ -0,0 +1,28 @@ +#include +#include +using namespace Rcpp; + + +// [[Rcpp::export]] +List mpolyMult(const List& e1, const List& e2) { + List res(e1.size() * e2.size()); + int c = 0; + for (int i = 0; i < e1.size(); i++) { + NumericVector v = e1[i]; + for (int j = 0; j < e2.size(); j++){ + NumericVector z = e2[j]; + res[c] = vcat(v, z); + c++; + } + } + return res; +} + +// [[Rcpp::export]] +List mpolyPow(const List& e1, int e2) { + List res = List::create(NumericVector::create(Named("coef", 1))); + for (int i = 0; i < e2; i++) { + res = mpolyMult(res, e1); + } + return res; +} From c47019a589e6035f7ea5ac50d8715fb7ba17615b Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:02:49 +0100 Subject: [PATCH 10/14] Make mpolyArithmetic.R use the Rcpp implementations --- R/mpolyArithmetic.R | 42 ++++++++++++++---------------------------- 1 file changed, 14 insertions(+), 28 deletions(-) diff --git a/R/mpolyArithmetic.R b/R/mpolyArithmetic.R index dacc6c0..acba3a8 100644 --- a/R/mpolyArithmetic.R +++ b/R/mpolyArithmetic.R @@ -94,36 +94,24 @@ NULL #' @rdname mpolyArithmetic #' @export `*.mpoly` <- function(e1, e2){ - ## allow for multiplication by a constant if(is.numeric(e1) && length(e1) == 1) e1 <- mpoly(list(c(coef = e1))) if(is.numeric(e2) && length(e2) == 1) e2 <- mpoly(list(c(coef = e2))) - - - + ## argument check stopifnot(is.mpoly(e1)) stopifnot(is.mpoly(e2)) - - - + ## multiply - list <- lapply(e1, function(v){ - lapply(e2, function(z){ - c(v, z) - }) - }) - - - + list <- mpolyMult(e1, e2) ## return - mpoly( unlist(list, recursive = FALSE) ) + mpoly( list ) } -l <- list( - c(x = 1, coef = 1, x = 1, coef = 1) -) +# l <- list( +# c(x = 1, coef = 1, x = 1, coef = 1) +# ) @@ -136,16 +124,14 @@ l <- list( #' @rdname mpolyArithmetic #' @export `^.mpoly` <- function(e1, e2){ - + if(!is.mpoly(e1)) stop('e1 must be of class mpoly.', call. = FALSE) - + if(!is.wholenumber(e2) || e2 < 0) stop('exponent must be a positive integer.', call. = FALSE) - - out <- mpoly(list(c(coef = 1))) - + if(e2 == 0) return(out) - - for(k in 1:e2) out <- out * e1 - - out + + out <- mpolyPow(e1, e2) + + mpoly(out) } From df7426f4062b1f6120a8845cfd2596be362e97ab Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:03:13 +0100 Subject: [PATCH 11/14] Add Rcpp dependency in NAMESPACE (generated by roxygen2) --- NAMESPACE | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 970f287..fe91ec7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -67,6 +67,8 @@ export(swap) export(totaldeg) export(tuples) export(vars) +import(Rcpp) +importFrom(Rcpp,sourceCpp) importFrom(ggplot2,aes) importFrom(ggplot2,ggplot) importFrom(orthopolynom,chebyshev.c.polynomials) @@ -101,3 +103,4 @@ importFrom(stringr,str_split) importFrom(stringr,str_sub) importFrom(stringr,str_trim) importFrom(tidyr,gather) +useDynLib(mpoly) From df2bb00673a68ccfef6121197c65af2901a95b77 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:03:30 +0100 Subject: [PATCH 12/14] Add Rcpp exports --- R/RcppExports.R | 35 +++++++++++++ src/RcppExports.cpp | 117 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 R/RcppExports.R create mode 100644 src/RcppExports.cpp diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..cf02a60 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,35 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +sortVars <- function(term, vars) { + .Call('_mpoly_sortVars', PACKAGE = 'mpoly', term, vars) +} + +getMonomialStrings <- function(poly) { + .Call('_mpoly_getMonomialStrings', PACKAGE = 'mpoly', poly) +} + +combineMonomials <- function(poly) { + .Call('_mpoly_combineMonomials', PACKAGE = 'mpoly', poly) +} + +tidyTerms <- function(poly, vars) { + .Call('_mpoly_tidyTerms', PACKAGE = 'mpoly', poly, vars) +} + +removeZeros <- function(poly) { + .Call('_mpoly_removeZeros', PACKAGE = 'mpoly', poly) +} + +reducePoly <- function(poly, vars) { + .Call('_mpoly_reducePoly', PACKAGE = 'mpoly', poly, vars) +} + +mpolyMult <- function(e1, e2) { + .Call('_mpoly_mpolyMult', PACKAGE = 'mpoly', e1, e2) +} + +mpolyPow <- function(e1, e2) { + .Call('_mpoly_mpolyPow', PACKAGE = 'mpoly', e1, e2) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..5d680e9 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,117 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace Rcpp; + +// sortVars +NumericVector sortVars(const NumericVector& term, const CharacterVector& vars); +RcppExport SEXP _mpoly_sortVars(SEXP termSEXP, SEXP varsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const NumericVector& >::type term(termSEXP); + Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); + rcpp_result_gen = Rcpp::wrap(sortVars(term, vars)); + return rcpp_result_gen; +END_RCPP +} +// getMonomialStrings +CharacterVector getMonomialStrings(const List& poly); +RcppExport SEXP _mpoly_getMonomialStrings(SEXP polySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); + rcpp_result_gen = Rcpp::wrap(getMonomialStrings(poly)); + return rcpp_result_gen; +END_RCPP +} +// combineMonomials +List combineMonomials(const List& poly); +RcppExport SEXP _mpoly_combineMonomials(SEXP polySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); + rcpp_result_gen = Rcpp::wrap(combineMonomials(poly)); + return rcpp_result_gen; +END_RCPP +} +// tidyTerms +List tidyTerms(const List& poly, const CharacterVector& vars); +RcppExport SEXP _mpoly_tidyTerms(SEXP polySEXP, SEXP varsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); + Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); + rcpp_result_gen = Rcpp::wrap(tidyTerms(poly, vars)); + return rcpp_result_gen; +END_RCPP +} +// removeZeros +List removeZeros(const List& poly); +RcppExport SEXP _mpoly_removeZeros(SEXP polySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); + rcpp_result_gen = Rcpp::wrap(removeZeros(poly)); + return rcpp_result_gen; +END_RCPP +} +// reducePoly +List reducePoly(const List& poly, const CharacterVector& vars); +RcppExport SEXP _mpoly_reducePoly(SEXP polySEXP, SEXP varsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); + Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); + rcpp_result_gen = Rcpp::wrap(reducePoly(poly, vars)); + return rcpp_result_gen; +END_RCPP +} +// mpolyMult +List mpolyMult(const List& e1, const List& e2); +RcppExport SEXP _mpoly_mpolyMult(SEXP e1SEXP, SEXP e2SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type e1(e1SEXP); + Rcpp::traits::input_parameter< const List& >::type e2(e2SEXP); + rcpp_result_gen = Rcpp::wrap(mpolyMult(e1, e2)); + return rcpp_result_gen; +END_RCPP +} +// mpolyPow +List mpolyPow(const List& e1, int e2); +RcppExport SEXP _mpoly_mpolyPow(SEXP e1SEXP, SEXP e2SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type e1(e1SEXP); + Rcpp::traits::input_parameter< int >::type e2(e2SEXP); + rcpp_result_gen = Rcpp::wrap(mpolyPow(e1, e2)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_mpoly_sortVars", (DL_FUNC) &_mpoly_sortVars, 2}, + {"_mpoly_getMonomialStrings", (DL_FUNC) &_mpoly_getMonomialStrings, 1}, + {"_mpoly_combineMonomials", (DL_FUNC) &_mpoly_combineMonomials, 1}, + {"_mpoly_tidyTerms", (DL_FUNC) &_mpoly_tidyTerms, 2}, + {"_mpoly_removeZeros", (DL_FUNC) &_mpoly_removeZeros, 1}, + {"_mpoly_reducePoly", (DL_FUNC) &_mpoly_reducePoly, 2}, + {"_mpoly_mpolyMult", (DL_FUNC) &_mpoly_mpolyMult, 2}, + {"_mpoly_mpolyPow", (DL_FUNC) &_mpoly_mpolyPow, 2}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_mpoly(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} From 1ef9efe6507f1049a68a3aaa2b8e5eb3e3a13246 Mon Sep 17 00:00:00 2001 From: Blaza Date: Sat, 30 Dec 2017 22:07:16 +0100 Subject: [PATCH 13/14] Remove excessive exports --- R/RcppExports.R | 20 --------------- src/RcppExports.cpp | 62 --------------------------------------------- src/mpoly.cpp | 7 ++--- 3 files changed, 2 insertions(+), 87 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index cf02a60..d351030 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,26 +1,6 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -sortVars <- function(term, vars) { - .Call('_mpoly_sortVars', PACKAGE = 'mpoly', term, vars) -} - -getMonomialStrings <- function(poly) { - .Call('_mpoly_getMonomialStrings', PACKAGE = 'mpoly', poly) -} - -combineMonomials <- function(poly) { - .Call('_mpoly_combineMonomials', PACKAGE = 'mpoly', poly) -} - -tidyTerms <- function(poly, vars) { - .Call('_mpoly_tidyTerms', PACKAGE = 'mpoly', poly, vars) -} - -removeZeros <- function(poly) { - .Call('_mpoly_removeZeros', PACKAGE = 'mpoly', poly) -} - reducePoly <- function(poly, vars) { .Call('_mpoly_reducePoly', PACKAGE = 'mpoly', poly, vars) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5d680e9..b796abb 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -5,63 +5,6 @@ using namespace Rcpp; -// sortVars -NumericVector sortVars(const NumericVector& term, const CharacterVector& vars); -RcppExport SEXP _mpoly_sortVars(SEXP termSEXP, SEXP varsSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const NumericVector& >::type term(termSEXP); - Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); - rcpp_result_gen = Rcpp::wrap(sortVars(term, vars)); - return rcpp_result_gen; -END_RCPP -} -// getMonomialStrings -CharacterVector getMonomialStrings(const List& poly); -RcppExport SEXP _mpoly_getMonomialStrings(SEXP polySEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); - rcpp_result_gen = Rcpp::wrap(getMonomialStrings(poly)); - return rcpp_result_gen; -END_RCPP -} -// combineMonomials -List combineMonomials(const List& poly); -RcppExport SEXP _mpoly_combineMonomials(SEXP polySEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); - rcpp_result_gen = Rcpp::wrap(combineMonomials(poly)); - return rcpp_result_gen; -END_RCPP -} -// tidyTerms -List tidyTerms(const List& poly, const CharacterVector& vars); -RcppExport SEXP _mpoly_tidyTerms(SEXP polySEXP, SEXP varsSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); - Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); - rcpp_result_gen = Rcpp::wrap(tidyTerms(poly, vars)); - return rcpp_result_gen; -END_RCPP -} -// removeZeros -List removeZeros(const List& poly); -RcppExport SEXP _mpoly_removeZeros(SEXP polySEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const List& >::type poly(polySEXP); - rcpp_result_gen = Rcpp::wrap(removeZeros(poly)); - return rcpp_result_gen; -END_RCPP -} // reducePoly List reducePoly(const List& poly, const CharacterVector& vars); RcppExport SEXP _mpoly_reducePoly(SEXP polySEXP, SEXP varsSEXP) { @@ -100,11 +43,6 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_mpoly_sortVars", (DL_FUNC) &_mpoly_sortVars, 2}, - {"_mpoly_getMonomialStrings", (DL_FUNC) &_mpoly_getMonomialStrings, 1}, - {"_mpoly_combineMonomials", (DL_FUNC) &_mpoly_combineMonomials, 1}, - {"_mpoly_tidyTerms", (DL_FUNC) &_mpoly_tidyTerms, 2}, - {"_mpoly_removeZeros", (DL_FUNC) &_mpoly_removeZeros, 1}, {"_mpoly_reducePoly", (DL_FUNC) &_mpoly_reducePoly, 2}, {"_mpoly_mpolyMult", (DL_FUNC) &_mpoly_mpolyMult, 2}, {"_mpoly_mpolyPow", (DL_FUNC) &_mpoly_mpolyPow, 2}, diff --git a/src/mpoly.cpp b/src/mpoly.cpp index 084d2d3..9ef1d54 100644 --- a/src/mpoly.cpp +++ b/src/mpoly.cpp @@ -52,7 +52,6 @@ struct compareVarOrder { }; -// [[Rcpp::export]] NumericVector sortVars(const NumericVector& term, const CharacterVector& vars) { std::vector names = term.names(); std::sort(names.begin(), names.end(), compareVarOrder(vars)); @@ -68,7 +67,6 @@ NumericVector sortVars(const NumericVector& term, const CharacterVector& vars) { } -// [[Rcpp::export]] CharacterVector getMonomialStrings(const List& poly) { CharacterVector res(poly.size()); @@ -94,7 +92,6 @@ CharacterVector getMonomialStrings(const List& poly) { } -// [[Rcpp::export]] List combineMonomials(const List& poly) { CharacterVector monomials = getMonomialStrings(poly); CharacterVector uniqueMonomials = unique(monomials); @@ -139,7 +136,7 @@ bool isSubset(CharacterVector sub, CharacterVector sup) { return true; } -// [[Rcpp::export]] + List tidyTerms(const List& poly, const CharacterVector& vars) { List res(poly.size()); @@ -166,7 +163,7 @@ List tidyTerms(const List& poly, const CharacterVector& vars) { return res; } -// [[Rcpp::export]] + List removeZeros(const List& poly){ // to avoid expanding the list, we'll first check how many zeros there are // and then proceed to populate the resulting list without zero elements From be9034fed8ac32c67f65d28304b009cc7eae9b06 Mon Sep 17 00:00:00 2001 From: Blaza Date: Fri, 6 Apr 2018 21:42:38 +0200 Subject: [PATCH 14/14] Add poly reducing in ^.mpoly This fixes the terrible performance of raising polynomials to a power. --- R/RcppExports.R | 4 ++-- R/mpolyArithmetic.R | 2 +- inst/include/helpers.h | 4 ++-- inst/include/reducepoly.h | 8 ++++++++ src/RcppExports.cpp | 9 +++++---- src/mpolyArithmetic.cpp | 5 +++-- 6 files changed, 21 insertions(+), 11 deletions(-) create mode 100644 inst/include/reducepoly.h diff --git a/R/RcppExports.R b/R/RcppExports.R index d351030..8e81456 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,7 +9,7 @@ mpolyMult <- function(e1, e2) { .Call('_mpoly_mpolyMult', PACKAGE = 'mpoly', e1, e2) } -mpolyPow <- function(e1, e2) { - .Call('_mpoly_mpolyPow', PACKAGE = 'mpoly', e1, e2) +mpolyPow <- function(e1, e2, vars) { + .Call('_mpoly_mpolyPow', PACKAGE = 'mpoly', e1, e2, vars) } diff --git a/R/mpolyArithmetic.R b/R/mpolyArithmetic.R index acba3a8..8eebeb7 100644 --- a/R/mpolyArithmetic.R +++ b/R/mpolyArithmetic.R @@ -131,7 +131,7 @@ NULL if(e2 == 0) return(out) - out <- mpolyPow(e1, e2) + out <- mpolyPow(e1, e2, vars(e1)) mpoly(out) } diff --git a/inst/include/helpers.h b/inst/include/helpers.h index 93b4e68..03e9f4f 100644 --- a/inst/include/helpers.h +++ b/inst/include/helpers.h @@ -1,5 +1,5 @@ -#ifndef mpoly_ADD_H -#define mpoly_ADD_H +#ifndef mpoly_HELPERS_H +#define mpoly_HELPERS_H using namespace Rcpp; diff --git a/inst/include/reducepoly.h b/inst/include/reducepoly.h new file mode 100644 index 0000000..f5124f0 --- /dev/null +++ b/inst/include/reducepoly.h @@ -0,0 +1,8 @@ +#ifndef mpoly_REDUCEPOLY_H +#define mpoly_REDUCEPOLY_H + +using namespace Rcpp; + +List reducePoly(const List& poly, const CharacterVector& vars); + +#endif diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b796abb..8856966 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -30,14 +30,15 @@ BEGIN_RCPP END_RCPP } // mpolyPow -List mpolyPow(const List& e1, int e2); -RcppExport SEXP _mpoly_mpolyPow(SEXP e1SEXP, SEXP e2SEXP) { +List mpolyPow(const List& e1, int e2, const CharacterVector& vars); +RcppExport SEXP _mpoly_mpolyPow(SEXP e1SEXP, SEXP e2SEXP, SEXP varsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const List& >::type e1(e1SEXP); Rcpp::traits::input_parameter< int >::type e2(e2SEXP); - rcpp_result_gen = Rcpp::wrap(mpolyPow(e1, e2)); + Rcpp::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); + rcpp_result_gen = Rcpp::wrap(mpolyPow(e1, e2, vars)); return rcpp_result_gen; END_RCPP } @@ -45,7 +46,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_mpoly_reducePoly", (DL_FUNC) &_mpoly_reducePoly, 2}, {"_mpoly_mpolyMult", (DL_FUNC) &_mpoly_mpolyMult, 2}, - {"_mpoly_mpolyPow", (DL_FUNC) &_mpoly_mpolyPow, 2}, + {"_mpoly_mpolyPow", (DL_FUNC) &_mpoly_mpolyPow, 3}, {NULL, NULL, 0} }; diff --git a/src/mpolyArithmetic.cpp b/src/mpolyArithmetic.cpp index 21cd445..7f32977 100644 --- a/src/mpolyArithmetic.cpp +++ b/src/mpolyArithmetic.cpp @@ -1,4 +1,5 @@ #include +#include #include using namespace Rcpp; @@ -19,10 +20,10 @@ List mpolyMult(const List& e1, const List& e2) { } // [[Rcpp::export]] -List mpolyPow(const List& e1, int e2) { +List mpolyPow(const List& e1, int e2, const CharacterVector& vars) { List res = List::create(NumericVector::create(Named("coef", 1))); for (int i = 0; i < e2; i++) { - res = mpolyMult(res, e1); + res = reducePoly(mpolyMult(res, e1), vars); } return res; }