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 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) diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..8e81456 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,15 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +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, vars) { + .Call('_mpoly_mpolyPow', PACKAGE = 'mpoly', e1, e2, vars) +} + 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 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) - - - diff --git a/R/mpolyArithmetic.R b/R/mpolyArithmetic.R index 9b4ce29..8eebeb7 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 } @@ -83,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) +# ) @@ -125,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, vars(e1)) + + mpoly(out) } 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 +} diff --git a/inst/include/helpers.h b/inst/include/helpers.h new file mode 100644 index 0000000..03e9f4f --- /dev/null +++ b/inst/include/helpers.h @@ -0,0 +1,24 @@ +#ifndef mpoly_HELPERS_H +#define mpoly_HELPERS_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/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/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 diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..8856966 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,56 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace 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, 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::traits::input_parameter< const CharacterVector& >::type vars(varsSEXP); + rcpp_result_gen = Rcpp::wrap(mpolyPow(e1, e2, vars)); + return rcpp_result_gen; +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, 3}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_mpoly(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} 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; +} diff --git a/src/mpoly.cpp b/src/mpoly.cpp new file mode 100644 index 0000000..9ef1d54 --- /dev/null +++ b/src/mpoly.cpp @@ -0,0 +1,214 @@ +#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; +}; + + +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; +} + + +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; +} + + +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; +} + + +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; +} + + +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; +} diff --git a/src/mpolyArithmetic.cpp b/src/mpolyArithmetic.cpp new file mode 100644 index 0000000..7f32977 --- /dev/null +++ b/src/mpolyArithmetic.cpp @@ -0,0 +1,29 @@ +#include +#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, const CharacterVector& vars) { + List res = List::create(NumericVector::create(Named("coef", 1))); + for (int i = 0; i < e2; i++) { + res = reducePoly(mpolyMult(res, e1), vars); + } + return res; +}