diff --git a/.development/check_as_cran.sh b/.development/check_as_cran.sh index 1b3b3ff..9ef78b5 100755 --- a/.development/check_as_cran.sh +++ b/.development/check_as_cran.sh @@ -1,5 +1,5 @@ #!/bin/bash -R CMD build ~/Documents/ast2ast -mv ./ast2ast_*.tar.gz ~/Documents/ast2ast/.development/ +R CMD build /home/konrad/Documents/GitHub/RProjects/ast2ast_supplement/ast2ast +mv ./ast2ast_*.tar.gz /home/konrad/Documents/GitHub/RProjects/ast2ast_supplement/ast2ast/.development R CMD check ./ast2ast_*.tar.gz --as-cran diff --git a/.development/install.sh b/.development/install.sh index 02e7391..99eceee 100755 --- a/.development/install.sh +++ b/.development/install.sh @@ -2,7 +2,9 @@ Rscript install.R -Rscript tracebackTest.R +Rscript testJacobian.R + +# Rscript tracebackTest.R # Rscript diff.R diff --git a/.development/testJacobian.R b/.development/testJacobian.R new file mode 100644 index 0000000..6e85f3c --- /dev/null +++ b/.development/testJacobian.R @@ -0,0 +1,397 @@ +f <- function(inp) { + if (inp == 0) { + b::double <- 4 + a <- c(1, 2, 3) + for (i in 1:length(a)) { + set_indep(a[i]) + a <- a * b + unset_indep(a[i]) + } + c <- get_deriv(a) + } else if (inp == 1) { + e::double <- 4 + for (i in 1:length(a)) { + set_indep(a[i]) + a <- a * b + unset_indep(a[i]) + } + c <- get_deriv(a) + } + return(c) +} +ast2ast::translate(f, verbose = TRUE) +stop() + + + +# Jacobian +# df1/dx1 ... df1/dxn +# dfm/dx1 ... dfm/dxn + +# https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant +f <- function(y, x, type_test) { + scal_d::double <- 0 # FIX: variable type declaration does not work in if + if (type_test == 1) { + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- (x[1]^2) * x[2] + y[2] <- 5 * x[1] + sin(x[2]) + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 2) { + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- x[1] * cos(x[2]) + y[2] <- x[1] * sin(x[2]) + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 3) { + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- x[1] * sin(x[2]) * cos(x[3]) + y[2] <- x[1] * sin(x[2]) * sin(x[3]) + y[3] <- x[1] * cos(x[2]) + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 4) { + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- x[1] + y[2] <- 5 * x[3] + y[3] <- 4 * x[2]^2 - 2 * x[3] + y[4] <- x[3] * sin(x[1]) + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 5) { + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- 5 * x[2] + y[2] <- 4 * x[1]^2 - 2 * sin(x[2] * x[3]) + y[3] <- x[2] * x[3] + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } + if (type_test == 6) { # Base case + y <- vector(length = 2) + # jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- x[1]^2 + sin(4) + y[2] <- x[2] * 7 + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 7) { # replacement at left side + y <- vector(length = 2) + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y1 <- y[1] + y1 <- x[1]^2 + y[2] <- x[2] * 7 + unset_indep(x[i]) + jac[, i] <- c(get_deriv(y1), get_deriv(y)[2]) + } + return(jac) + } else if (type_test == 8) { # replacement at right side + y <- rep(0, 2) + a <- 8 + jac <- matrix(0, length(y), length(x)) + for (i in 1:length(x)) { + set_indep(x[i]) + y1 <- y[1] + y1 <- x[1]^2 + y[2] <- x[2] * a + unset_indep(x[i]) + jac[, i] <- c(get_deriv(y1), get_deriv(y[2])) + } + return(jac) + } else if (type_test == 9) { # replacement at rs + y <- vector(length = length(x)) + jac <- matrix(0, length(y), length(x)) + + for (i in 1:length(x)) { + set_indep(x[i]) + a <- 4 * x[1] + y[1] <- (x[1]^2) * a + y[2] <- x[2] * a + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + + return(jac) + } else if (type_test == 10) { # replacement at rs & ls + y <- vector(length = length(x)) + jac <- matrix(0, length(y), length(x)) + + for (i in 1:length(x)) { + set_indep(x[i]) + a <- 4 * x[1] + y1 <- y[1] + y2 <- y[2] + y1 <- (x[1]^2) * a + y2 <- x[2] * a + unset_indep(x[i]) + jac[, i] <- c(get_deriv(y1), get_deriv(y2)) + } + + return(jac) + } else if (type_test == 11) { # if and replace at ls + + y <- vector(length = length(x)) + jac <- matrix(0, length(y), length(x)) + t <- 3.5 + for (i in 1:length(x)) { + set_indep(x[i]) + y[1] <- 2 * (x[1]^3) # 6x^2 = 24, 0 + if (t < 3) { + y[2] <- x[2]^2 + } else if (t < 5) { + y[2] <- x[2]^3 # 3x^2 = 0, 75 + } else { + y[2] <- x[2]^4 + } + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } else if (type_test == 12) { # if and replace at ls & rs + y <- vector(length = length(x)) + jac <- matrix(0, length(y), length(x)) + t::double <- 2.5 + for (i in 1:length(x)) { + set_indep(x[i]) + y2 <- y[2] + a <- x[1] * 3 + y[1] <- 2 * x[1]^3 + if (t < 3) { + y2 <- x[2]^2 + } else if (t < 5) { + y2 <- x[2]^3 + } else { + y2 <- x[2]^4 * a + } + unset_indep(x[i]) + jac[, i] <- c(get_deriv(y[1]), get_deriv(y2)) + } + return(jac) + } else if (type_test == 13) { + # function from wikipedia + # https://de.wikipedia.org/wiki/Jacobi-Matrix + b <- vector(length = 2) + jac <- matrix(0, 2, 3) + a <- x + + for (i in 1:length(a)) { + set_indep(a[i]) + b[1] <- a[1]^2 + a[2]^2 + a[3] * sin(a[1]) + b[2] <- a[3]^2 + a[3] * sin(a[2]) + unset_indep(a[i]) + jac[, i] <- get_deriv(b) + } + + return(jac) + } else if (type_test == 14) { # test cmr + t::double <- 4 + times <- x + values <- y + x <- c(2, 5) + y <- vector(length = 2) + # FIX: vector(length(x)) this results in compile error should be catched by function mapping + jac <- matrix(0, length(x), length(y)) + for (i in 1:length(x)) { + set_indep(x[i]) + a <- cmr(t, times, values) + x1 <- x[1] + x2 <- x[2] + y1 <- y[1] + y2 <- y[2] + y1 <- a * x1 + y2 <- a + x2 + unset_indep(x[i]) + jac[, i] <- c(get_deriv(y1), get_deriv(y2)) + } + return(jac) + } else if (type_test == 15) { # test scalar values + y <- rep(0, length(x)) + jac <- matrix(0, length(x), length(y)) + for (i in 1:length(x)) { + set_indep(x[i]) + scal_d <- x[1] + y[1] <- scal_d * scal_d + y[2] <- x[1] * x[1] + unset_indep(x[i]) + jac[, i] <- get_deriv(y) + } + return(jac) + } + return() +} + +jac1 <- function(x) { + jac <- matrix(0, 2, 2) + jac[1, ] <- c(2 * x[1] * x[2], x[1]^2) + jac[2, ] <- c(5, cos(x[2])) + return(jac) +} + +jac2 <- function(x) { + jac <- matrix(0, 2, 2) + jac[1, ] <- c(cos(x[2]), -x[1] * sin(x[2])) + jac[2, ] <- c(sin(x[2]), x[1] * cos(x[2])) + return(jac) +} + +jac3 <- function(x) { + jac <- matrix(0, 3, 3) + jac[1, ] <- c( + sin(x[2]) * cos(x[3]), + x[1] * cos(x[2]) * cos(x[3]), + -x[1] * sin(x[2]) * sin(x[3]) + ) + jac[2, ] <- c( + sin(x[2]) * sin(x[3]), + x[1] * cos(x[2]) * sin(x[3]), + x[1] * sin(x[2]) * cos(x[3]) + ) + jac[3, ] <- c(cos(x[2]), -x[1] * sin(x[2]), 0) + return(jac) +} + +jac4 <- function(x) { + jac <- matrix(0, 4, 3) + jac[1, ] <- c(1, 0, 0) + jac[2, ] <- c(0, 0, 5) + jac[3, ] <- c(0, 8 * x[2], -2) + jac[4, ] <- c(x[3] * cos(x[1]), 0, sin(x[1])) + return(jac) +} + +jac5 <- function(x) { + jac <- matrix(0, 3, 3) + jac[1, ] <- c(0, 5, 0) + jac[2, ] <- c(8 * x[1], -2 * x[3] * cos(x[2] * x[3]), -2 * x[2] * cos(x[2] * x[3])) + jac[3, ] <- c(0, x[3], x[2]) + return(jac) +} + +fcpp <- ast2ast::translate(f, + types_of_args = c("double", "double", "int"), + data_structures = c("vector", "vector", "scalar"), + handle_inputs = c("copy", "copy", "copy"), + verbose = TRUE +) + +x <- c(1, 2) +y <- numeric(2) +expect_equal(fcpp(y, x, 1L), jac1(x)) + +x <- c(1, 2) +y <- numeric(2) +expect_equal(fcpp(y, x, 2L), jac2(x)) + +x <- c(1, 2, 3) +y <- numeric(3) +expect_equal(fcpp(y, x, 3L), jac3(x)) + +x <- c(1, 2, 3) +y <- numeric(4) +expect_equal(fcpp(y, x, 4L), jac4(x)) + +x <- c(1, 2, 3) +y <- numeric(3) +expect_equal(fcpp(y, x, 5L), jac5(x)) + + +# basic case +d <- c(4, 0, 0, 7) +nrow <- 2 +ncol <- 2 +x <- c(2, 5) +y <- numeric(2) +expect_equal(fcpp(y, x, 6L), matrix(d, nrow, ncol)) + +# replacement at left side +x <- c(2, 5) +res <- fcpp(y, x, 7L) +d <- c(4, 0, 0, 7) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# replacement at right side +x <- c(2, 5) +res <- fcpp(y, x, 8L) +d <- c(4, 0, 0, 8) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# replacement at rs +x <- c(2, 5) +res <- fcpp(y, x, 9L) +d <- c(48, 20, 0, 8) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# replacement at ls & rs +x <- c(2, 5) +res <- fcpp(y, x, 10L) +d <- c(48, 20, 0, 8) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# replacement at ls & rs +x <- c(2, 5) +res <- fcpp(y, x, 11L) +d <- c(24, 0, 0, 75) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# if & replacement at ls & rs +x <- c(2, 5) +res <- fcpp(y, x, 12L) +d <- c(24, 0, 0, 10) +nrow <- 2 +ncol <- 2 +expect_equal(res, matrix(d, nrow, ncol)) + +# wikipedia +# https://de.wikipedia.org/wiki/Jacobi-Matrix +x <- c(pi, pi, pi) +y <- numeric(2) +res <- fcpp(y, x, 13L) +d <- c(pi, 0, pi * 2, -pi, 1.224647e-16, pi * 2) +nrow <- 2 +ncol <- 3 +expect_equal(round(res, 2), round(matrix(d, nrow, ncol), 2)) + +# cmr +times <- seq(0, 24, 2) +values <- runif(13) +res <- fcpp(values, times, 14L) +expect_equal(res, matrix(c(values[3], 0, 0, 1), 2, 2)) + +# scalar values +x <- c(1, 2) +res <- fcpp(y, x, 15L) +expect_equal(res, matrix(c(0, 2, 0, 0), 2, 2)) diff --git a/R/AstClass.R b/R/AstClass.R index 1ec63cf..981881a 100644 --- a/R/AstClass.R +++ b/R/AstClass.R @@ -29,8 +29,6 @@ astClass <- R6::R6Class("astClass", }, error = function(err) { e <- err - # print_traceback() - # print_debug_info() e } ) diff --git a/R/CodelineClass.R b/R/CodelineClass.R index a9aa361..d39c24a 100644 --- a/R/CodelineClass.R +++ b/R/CodelineClass.R @@ -53,7 +53,7 @@ generic_fcts <- function() { c( "+", "-", "*", "/", "if", "else if", "else", "{", "(", - "==", "!=", ">", ">=", "<", "<=", "vector", + "==", "!=", ">", ">=", "<", "<=", "rep", "::", "matrix", "length", "dim", "cmr", "exp", "at", "&&", "||", "Rf_ScalarReal", "cpp2R", @@ -134,6 +134,10 @@ LC <- R6::R6Class("LC", } else if (length(sexp) == 1) { self$return_variable <- NULL } + } else if (as.name("vector") == fct) { + self$check_assign_subset <- FALSE + p <- vector_class$new(sexp, self$namespace_etr) + sexp <- p$convert(self$PF) } else if (as.name("[") == fct) { p <- subset$new(sexp, self$check_assign_subset, self$namespace_etr) sexp <- p$convert(self$PF) diff --git a/R/Compiling.R b/R/Compiling.R index 017c69d..c6c97ea 100644 --- a/R/Compiling.R +++ b/R/Compiling.R @@ -40,7 +40,7 @@ compiler_a2a <- function(fct_code, R_fct, attributes(fct_ret) <- list(class = "XPtr") }, error = function(e) { - print("Sorry compilation failed!") + color_print(43, "Sorry compilation failed!") } ) } else { @@ -58,14 +58,14 @@ compiler_a2a <- function(fct_code, R_fct, fct_ret <- env[[name_f]] }, error = function(e) { - print("Sorry compilation failed!") + color_print(43, "Sorry compilation failed!") } ) Sys.unsetenv("PKG_CXXFLAGS") if (verbose == TRUE) { - cat(fct) + print_with_line_numbers(fct) } } diff --git a/R/Differentiate.R b/R/Differentiate.R index 756c161..60f4493 100755 --- a/R/Differentiate.R +++ b/R/Differentiate.R @@ -160,7 +160,10 @@ diff_built_in_function_call <- lift(function(expr, fl) { od <- list() counter <- 1 for (i in seq_along(args)) { - od[[counter]] <- pryr::modify_lang(outer_deriv, substi, bquote(.(deriv_args[[i]])), args[[i]]) + od[[counter]] <- pryr::modify_lang( + outer_deriv, substi, + bquote(.(deriv_args[[i]])), args[[i]] + ) counter <- counter + 1 } outer_deriv <- od @@ -212,6 +215,9 @@ diff_call <- lift(function(expr, fl) { arg1 <- call_arg(expr, 1) arg2 <- call_arg(expr, 2) e <- call_name(expr) + if (deparse(e) == "etr::i2d") { + return(0) + } if (is.name(e)) { if (e == "+") { diff_addition(expr, fl) diff --git a/R/FctDerivPairs.R b/R/FctDerivPairs.R index 5c4dcf0..afc8d9c 100755 --- a/R/FctDerivPairs.R +++ b/R/FctDerivPairs.R @@ -167,6 +167,15 @@ numeric_deriv <- function(x) numeric(x) rep_deriv <- function(x, y) rep(x, y) matrix_deriv <- function(val = 0, x, y) matrix(val, x, y) +cmr_deriv <- function(a, b, c, d) cmr(a, b, c, d) +get_deriv_deriv <- function(a) get_deriv(a) +colon_deriv <- function(a, b) a:b +length_deriv <- function(a) length(a) +dim_deriv <- function(a) dim(a) +set_indep_deriv <- function(a) set_indep(a) +unset_indep_deriv <- function(a) unset_indep(a) +i2d_deriv <- function(a) etr::i2d(a) + init_fct_list <- function() { @@ -186,13 +195,14 @@ init_fct_list <- function() { f <- add_fct(f, "vector", vector_deriv, TRUE) f <- add_fct(f, "rep", rep_deriv, TRUE) f <- add_fct(f, "matrix", matrix_deriv, TRUE) - f <- add_fct(f, "c", matrix_deriv, TRUE) - f <- add_fct(f, "cmr", matrix_deriv, TRUE) - f <- add_fct(f, "get_deriv", matrix_deriv, TRUE) - f <- add_fct(f, ":", matrix_deriv, TRUE) - f <- add_fct(f, "length", matrix_deriv, TRUE) - f <- add_fct(f, "dim", matrix_deriv, TRUE) - f <- add_fct(f, "set_indep", matrix_deriv, TRUE) - f <- add_fct(f, "unset_indep", matrix_deriv, TRUE) + f <- add_fct(f, "c", c_deriv, TRUE) + f <- add_fct(f, "cmr", cmr_deriv, TRUE) + f <- add_fct(f, "get_deriv", get_deriv_deriv, TRUE) + f <- add_fct(f, ":", colon_deriv, TRUE) + f <- add_fct(f, "length", length_deriv, TRUE) + f <- add_fct(f, "dim", dim_deriv, TRUE) + f <- add_fct(f, "set_indep", set_indep_deriv, TRUE) + f <- add_fct(f, "unset_indep", unset_indep_deriv, TRUE) + f <- add_fct(f, "etr::i2d", i2d_deriv, TRUE) return(f) } diff --git a/R/NodeClasses.R b/R/NodeClasses.R index b1a2f7e..68a5702 100644 --- a/R/NodeClasses.R +++ b/R/NodeClasses.R @@ -216,7 +216,7 @@ generic <- R6::R6Class("generic", inherit = PC, public = list( change_code = function() { - self$replace_int() # check TRUE FALSE? + self$replace_int() }, convert = function(var) { self$replace_INF() @@ -238,6 +238,30 @@ generic <- R6::R6Class("generic", ) ) +#' @import R6 +vector_class <- R6::R6Class("vector_class", + inherit = PC, + public = list( + convert = function(var) { + self$replace_INF() + self$replace_NA() + self$replace_TF() + self$oaf(var) + + ret <- list() + if (deparse(self$name_fct) %in% self$namespace_etr) { + ret[[1]] <- as.name(paste0("etr::", self$name_fct)) + } else { + ret[[1]] <- self$name_fct + } + + ret <- c(ret, self$arguments) + return(ret) + } + ) +) + + #' @import R6 assign <- R6::R6Class("assign", inherit = PC, diff --git a/R/Utils.R b/R/Utils.R index 56767d4..240a7b6 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -22,6 +22,23 @@ color_print <- function(col, ...) { cat(paste0("\033[0;", col, "m", txt, "\033[0m", "\n")) } +add_line_numbers <- function(text) { + sprintf("%d: %s", seq_along(lines), lines) +} + +color_print_wrong_code <- function(col, ...) { + txt <- c(...) + txt <- add_line_numbers(txt) + txt <- paste(txt, collapse = " ") + cat(paste0("\033[0;", col, "m", txt, "\033[0m", "\n")) +} + +print_with_line_numbers <- function(text) { + lines <- strsplit(text, "\n")[[1]] + numbered_lines <- sprintf("%d: %s", seq_along(lines), lines) + cat(paste(numbered_lines, collapse = "\n")) +} + print_traceback <- function() { tb <- sys.calls() if (length(tb) > 0) {