Expanding a polynomial with ‘caracas’, part 2
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Last month, I posted an article showing a way to expand a polynomial in R when the coefficients of the polynomial contain some literal values, with the help of the caracas package. Today I wanted to apply it with a polynomial expression having about 500 characters, and highly factorized. The method took more than 30 minutes, so I looked for a more efficient one.
Thanks to some help on StackOverflow, I came to the following method which is more efficient. It consists in splitting the expression according to its additive terms and to work on each term, instead of expanding the whole polynomial. In the example below I take the polynomial expression defining the isosurface equation of a Solid Möbius strip.
library(caracas)
sympy <- get_sympy()
# define the variables x,y,z and the constants a,b
def_sym(x, y, z, a, b)
# define expression
expr <- sympy$parse_expr(
"((x*x+y*y+1)*(a*x*x+b*y*y)+z*z*(b*x*x+a*y*y)-2*(a-b)*x*y*z-a*b*(x*x+y*y))**2-4*(x*x+y*y)*(a*x*x+b*y*y-x*y*z*(a-b))**2"
)
# extraction of monomials in the 'povray' list
povray <- list()
terms <- sympy$Add$make_args(expr)
for(term in terms){
f <- term$expand()
fterms <- sympy$Add$make_args(f)
for(fterm in fterms){
decomp <- fterm$as_coeff_mul(x$pyobj, y$pyobj, z$pyobj)
coef <- decomp[[1]]
mono <- decomp[[2]]
polexpr <- sympy$Mul$fromiter(mono)
poly <- polexpr$as_poly(x$pyobj, y$pyobj, z$pyobj)
degree <- toString(poly$monoms()[[1]])
if(degree %in% names(povray)){
povray[[degree]] <- sympy$Add(povray[[degree]], coef)
}else{
povray[[degree]] <- coef
}
}
}
polynomial <- vapply(names(povray), function(degree){
coeff <- povray[[degree]] |>
gsub("([ab])\\*\\*(\\d+)", "pow(\\1,\\2)", x = _)
sprintf("xyz(%s): %s,", degree, coeff)
}, character(1L))
cat(polynomial, sep = "\n", file = "SolidMobiusStrip.txt")
At the last step I use gsub to replace the powers like
a**2 to their POV-Ray syntax pow(a,2). The
above code writes this POV-Ray code:
xyz(4, 0, 0): pow(a,2)*pow(b,2) - 2*pow(a,2)*b + pow(a,2), xyz(8, 0, 0): pow(a,2), xyz(0, 4, 0): pow(a,2)*pow(b,2) - 2*a*pow(b,2) + pow(b,2), xyz(0, 8, 0): pow(b,2), xyz(6, 0, 0): -2*pow(a,2)*b - 2*pow(a,2), xyz(0, 6, 0): -2*a*pow(b,2) - 2*pow(b,2), xyz(4, 4, 0): pow(a,2) + 4*a*b + pow(b,2), xyz(0, 4, 4): pow(a,2), xyz(4, 0, 4): pow(b,2), xyz(4, 2, 0): -4*pow(a,2)*b - 2*pow(a,2) - 2*a*pow(b,2) - 4*a*b, xyz(6, 2, 0): 2*pow(a,2) + 2*a*b, xyz(2, 4, 0): -2*pow(a,2)*b - 4*a*pow(b,2) - 4*a*b - 2*pow(b,2), xyz(2, 6, 0): 2*a*b + 2*pow(b,2), xyz(1, 3, 3): -4*pow(a,2) + 4*a*b, xyz(3, 1, 1): 4*pow(a,2)*b - 4*pow(a,2) - 4*a*pow(b,2) + 4*a*b, xyz(5, 1, 1): 4*pow(a,2) - 4*a*b, xyz(3, 3, 1): 4*pow(a,2) - 4*pow(b,2), xyz(2, 2, 0): 2*pow(a,2)*pow(b,2) - 2*pow(a,2)*b - 2*a*pow(b,2) + 2*a*b, xyz(4, 0, 2): -2*a*pow(b,2) + 2*a*b, xyz(0, 4, 2): -2*pow(a,2)*b + 2*a*b, xyz(6, 0, 2): 2*a*b, xyz(0, 6, 2): 2*a*b, xyz(2, 4, 2): -2*pow(a,2) + 10*a*b - 2*pow(b,2), xyz(4, 2, 2): -2*pow(a,2) + 10*a*b - 2*pow(b,2), xyz(1, 3, 1): 4*pow(a,2)*b - 4*a*pow(b,2) - 4*a*b + 4*pow(b,2), xyz(1, 5, 1): 4*a*b - 4*pow(b,2), xyz(3, 1, 3): -4*a*b + 4*pow(b,2), xyz(2, 2, 2): -2*pow(a,2)*b + 6*pow(a,2) - 2*a*pow(b,2) - 8*a*b + 6*pow(b,2), xyz(2, 2, 4): 2*a*b,
This is very fast for this example, but it still took 20 minutes with my case, which is a slight modification of an animation by ‘ICN5D’; here it is:
The difference with the original animation is that this one uses an isoclinic rotation for the animation.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.