R/poly_solve.R

In dkahle/bertini: Interface to 'Bertini'

#' Solve a System of Polynomial Equations
#'
#' [poly_solve()] solves a system of polynomial equations, specifiable in any
#' of several ways.
#'
#' @param lhs a mpolyList or character vector of left hand sides
#' @param rhs a mpolyList or character vector of right hand sides
#' @param varorder variable order (see examples)
#' @param ... stuff to pass to bertini
#' @return an object of class bertini
#' @export poly_solve
#' @seealso [variety()], [bertini()]
#' @examples
#'
#' if (has_bertini()) {
#'
#' # it can solve linear systems
#' # (here where the line y = x intersects y = 2 - x)
#' poly_solve(c("y", "y"), c("x", "2 - x"), c("x", "y"))
#'
#' # or nonlinear systems
#' poly_solve(c("y", "y"), c("x^2", "2 - x^2"), c("x", "y"))
#'
#' # perhaps an easier specification is equations themselves
#' # with either the " = " or " == " specifications
#' # varorder is used to order the solutions returned
#' poly_solve(c("y = x^2", "y = 2 - x^2"), varorder = c("x", "y"))
#' poly_solve(c("y == x^2", "y == 2 - x^2"), varorder = c("x", "y"))
#'
#'
#' # mpoly objects can be given instead of character strings
#' lhs <- mp(c("y - (2 - x)", "x y"))
#' rhs <- mp(c("0","0"))
#' poly_solve(lhs, rhs, varorder = c("x", "y"))
#'
#' # if no default right hand side is given, and no "=" or "==" is found,
#' # rhs is taken to be 0's.
#' # below is where the lines y = x and y = -x intersect the unit circle
#' poly_solve(c("(y - x) (y + x)", "x^2 + y^2 - 1"))
#'
#' # the output object is a bertini object
#' out <- poly_solve(c("(y - x) (y + x)", "x^2 + y^2 - 1"))
#' str(out, 1)
#'
#' # here is the code that was run :
#' cat(out$bertini_code)
#'
#' # the finite and real solutions:
#' out$finite_solutions
#' out$real_finite_solutions
#'
#'
#'
#' # (known priting issue)
#' # example from Riccomagno (2008), p. 399
#' # poly_solve(c(
#' #   "x (x - 2) (x - 4) (x - 3)",
#' #   "(y - 4) (y - 2) y",
#' #   "(y - 2) (x + y - 4)",
#' #   "(x - 3) (x + y - 4)"
#' # ))
#'
#' }
#'
poly_solve <- function(lhs, rhs, varorder, ...){

  # this should use mpoly::eq_mp in the next version of bertini
  # that function is in mpoly_1_1_1; mpoly_1_1_0 is currently (02/24/19) on
  # cran. the first release of bertini 0.0.1 will depend on mpoly_1_1_0

  if(missing(rhs)){

    if(is.character(lhs)){
      if(all(str_detect(lhs, "==?"))){
        split <- str_split(lhs, " ==? ")
        lhs <- vapply(split, function(x) x[1], character(1))
        rhs <- vapply(split, function(x) x[2], character(1))
      } else {
        rhs <- mp(rep("0", length(lhs)))
        if(length(rhs) == 1){
          rhs <- list(rhs)
          class(rhs) <- "mpolyList"
        }
      }
    } else {
      if(is.mpoly(lhs)){
        rhs <- list(mp("0"))
        class(rhs) <- "mpolyList"
      } else if(is.mpolyList(lhs)){
        rhs <- mp(rep("0", length(lhs)))
      }
    }
  }

  if(is.character(lhs)) lhs <- mp(lhs)

  if(is.mpoly(lhs)) lhs <- structure(list(lhs), class = "mpolyList")

  if(!missing(rhs)){
    if(is.character(rhs)) rhs <- mp(rhs)
    if(is.mpoly(rhs)) rhs <- structure(list(rhs), class = "mpolyList")
  }

  stopifnot(is.mpolyList(lhs) && is.mpolyList(rhs))

  gens <- lhs - rhs


  # sort out variables
  vars <- vars(gens)

  if(!missing(varorder) && !all(sort(vars) == sort(varorder))) stop(
    "If varorder is provided, it must contain all of the variables.",
    call. = FALSE
  )
  if(!missing(varorder)) vars <- varorder


  # format code
  fun_names <- str_c("f", 1:length(gens))
  funs <- print(gens, stars = TRUE, silent = TRUE)
  funs <- str_replace_all(funs, "  ", " ")
  funs <- str_replace_all(funs, "\\*\\*", "^")

  # glue code together
  code <- glue(
"INPUT

variable_group {paste(vars, collapse = ', ')};
function {paste(fun_names, collapse = ', ')};
{paste(fun_names, str_c(funs, ';'), sep = ' = ', collapse = '\n')}

END;"
  )

  bertini(code, ...)
}