From 70cd4c03209c4571799e335895fe8fb2c14d77f4 Mon Sep 17 00:00:00 2001 From: Ben Bolker Date: Thu, 2 May 2024 16:16:15 -0400 Subject: [PATCH 1/4] first pass at safe_power() engine fn --- R/engine_functions.R | 1 + src/macpan2.cpp | 10 ++++++++++ 2 files changed, 11 insertions(+) diff --git a/R/engine_functions.R b/R/engine_functions.R index 7d0259810..d0858686a 100644 --- a/R/engine_functions.R +++ b/R/engine_functions.R @@ -816,4 +816,5 @@ #' @aliases reulermultinom #' @aliases round #' @aliases pgamma +#' @aliases safe_power NULL diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 7060c5c0a..9bba8d0d1 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -112,6 +112,8 @@ enum macpan2_func , MP2_EULER_MULTINOM_SIM = 48 // fwrap,fail: reulermultinom(size, rate, delta_t) , MP2_ROUND = 49 // fwrap,null: round(x) , MP2_PGAMMA = 50 // fwrap,fail: pgamma(q, shape, scale) + // `^`(x, y), but with 0^0 defined as 0 + , MP2_SAFEPOWER = 51 // fwrap,fail: safe_power(x,y) }; enum macpan2_meth @@ -1288,6 +1290,14 @@ class ExprEvaluator return pow(args[0].array(), args[1].array()).matrix(); // return args[0].pow(args[1].coeff(0,0)); + case MP2_SAFEPOWER: // ^ + #ifdef MP_VERBOSE + std::cout << "I'm being lazy" << std::endl; + #endif + // ?? are the CondExpRel functions vectorized ... ??? + // CondExpRel(left, right, if_true, if_false) + return CppAD::CondExpLe(args[0].array(), Type(0), Type(0), pow(args[0].array(), args[1].array()).matrix()); + // #' ## Unary Elementwise Math // #' // #' ### Functions From 07fc4224061acdcc4e469199637768954438a080 Mon Sep 17 00:00:00 2001 From: Ben Bolker Date: Thu, 2 May 2024 17:13:06 -0400 Subject: [PATCH 2/4] add correct SAFEPOWER code --- man/engine_functions.Rd | 1 + src/macpan2.cpp | 38 +++++++++++++++++++++++++++++++------- 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/man/engine_functions.Rd b/man/engine_functions.Rd index e574d9a48..dc35d531e 100644 --- a/man/engine_functions.Rd +++ b/man/engine_functions.Rd @@ -52,6 +52,7 @@ \alias{reulermultinom} \alias{round} \alias{pgamma} +\alias{safe_power} \title{Engine Functions} \description{ Functions currently supported by the C++ TMB engine diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 9bba8d0d1..672ddfb50 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -1290,13 +1290,37 @@ class ExprEvaluator return pow(args[0].array(), args[1].array()).matrix(); // return args[0].pow(args[1].coeff(0,0)); - case MP2_SAFEPOWER: // ^ - #ifdef MP_VERBOSE - std::cout << "I'm being lazy" << std::endl; - #endif - // ?? are the CondExpRel functions vectorized ... ??? - // CondExpRel(left, right, if_true, if_false) - return CppAD::CondExpLe(args[0].array(), Type(0), Type(0), pow(args[0].array(), args[1].array()).matrix()); + case MP2_SAFEPOWER: // SAFE_POWER, equivalent to (ifelse(x==0, 0, x^y)) + + args = args.recycle_for_bin_op(); + err_code = args.get_error_code(); + switch (err_code) { + case 201: + SetError(err_code, "The two operands do not have the same number of columns", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + case 202: + SetError(err_code, "The two operands do not have the same number of rows", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + case 203: + SetError(err_code, "The two operands do not have the same number of columns or rows", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + } + m1 = args[0]; + m2 = args[1]; + rows = m1.rows(); + cols = m1.cols(); + m = matrix::Zero(rows, cols); + for (int i; i < rows; i++) { + for (int j; j < cols; j++) { + m.coeffRef(i, j) = CppAD::CondExpEq( + m1.coeff(i, j), + m2.coeff(i, j), + Type(0), + pow(m1.coeff(i, j), m2.coeff(i, j)) + ); + } + } + return m; // #' ## Unary Elementwise Math // #' From 6965bf4277e5d936cfee8e3edca163975e54424b Mon Sep 17 00:00:00 2001 From: Ben Bolker Date: Thu, 2 May 2024 17:34:31 -0400 Subject: [PATCH 3/4] update workflows to checkout@v4 --- .github/workflows/R-CMD-check.yaml | 2 +- .github/workflows/pkgdown.yaml | 2 +- .github/workflows/test-coverage.yaml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index cd9685142..bf5ea8fe2 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -30,7 +30,7 @@ jobs: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index d6ffea3a2..02c40ef26 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -18,7 +18,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 9b1b5898b..3f9d3a27b 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -13,7 +13,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 with: From fdfef676a7b403485ce1458e930d2ebdc88cbc1d Mon Sep 17 00:00:00 2001 From: stevencarlislewalker Date: Fri, 3 May 2024 05:03:19 -0400 Subject: [PATCH 4/4] first pass at #200 --- R/enum.R | 1 + misc/dev/dev.cpp | 36 ++++++++++++++++++++++++++++++++++++ src/macpan2.cpp | 16 +++++++++------- tests/testthat/test-binop.R | 29 +++++++++++++++++++++++++++++ 4 files changed, 75 insertions(+), 7 deletions(-) diff --git a/R/enum.R b/R/enum.R index 9b8183609..4bfa5bc6a 100644 --- a/R/enum.R +++ b/R/enum.R @@ -50,6 +50,7 @@ valid_func_sigs = c( , "fwrap,fail: reulermultinom(size, rate, delta_t)" , "fwrap,null: round(x)" , "fwrap,fail: pgamma(q, shape, scale)" + , "fwrap,fail: safe_power(x,y)" ) process_enum = function(x) { RE = "(null|fail|binop|fwrap|bwrap|pwrap)[ ]*,[ ]*(null|fail|binop|fwrap|bwrap|pwrap)[ ]*:[ ]*\\`?([^`]*)\\`?\\((.*)(\\,.*)*\\)" diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp index a63d07503..752250bc8 100644 --- a/misc/dev/dev.cpp +++ b/misc/dev/dev.cpp @@ -111,6 +111,7 @@ enum macpan2_func , MP2_EULER_MULTINOM_SIM = 48 // fwrap,fail: reulermultinom(size, rate, delta_t) , MP2_ROUND = 49 // fwrap,null: round(x) , MP2_PGAMMA = 50 // fwrap,fail: pgamma(q, shape, scale) + , MP2_SAFEPOWER = 51 // fwrap,fail: safe_power(x,y) }; enum macpan2_meth @@ -1286,6 +1287,41 @@ class ExprEvaluator #endif return pow(args[0].array(), args[1].array()).matrix(); // return args[0].pow(args[1].coeff(0,0)); + + case MP2_SAFEPOWER: // SAFE_POWER, equivalent to (ifelse(x==0, 0, x^y)) + if (n != 2) { + SetError(err_code, "safe_power requires exactly two arguments", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + } + args = args.recycle_for_bin_op(); + err_code = args.get_error_code(); + switch (err_code) { + case 201: + SetError(err_code, "The two operands do not have the same number of columns", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + case 202: + SetError(err_code, "The two operands do not have the same number of rows", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + case 203: + SetError(err_code, "The two operands do not have the same number of columns or rows", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + } + m1 = args[0]; + m2 = args[1]; + rows = m1.rows(); + cols = m1.cols(); + m = matrix::Zero(rows, cols); + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + m.coeffRef(i, j) = CppAD::CondExpEq( + m1.coeff(i, j), + Type(0), + Type(0), + pow(m1.coeff(i, j), m2.coeff(i, j)) + ); + } + } + return m; // #' ## Unary Elementwise Math // #' diff --git a/src/macpan2.cpp b/src/macpan2.cpp index 672ddfb50..8e13a2344 100644 --- a/src/macpan2.cpp +++ b/src/macpan2.cpp @@ -112,7 +112,6 @@ enum macpan2_func , MP2_EULER_MULTINOM_SIM = 48 // fwrap,fail: reulermultinom(size, rate, delta_t) , MP2_ROUND = 49 // fwrap,null: round(x) , MP2_PGAMMA = 50 // fwrap,fail: pgamma(q, shape, scale) - // `^`(x, y), but with 0^0 defined as 0 , MP2_SAFEPOWER = 51 // fwrap,fail: safe_power(x,y) }; @@ -1289,9 +1288,12 @@ class ExprEvaluator #endif return pow(args[0].array(), args[1].array()).matrix(); // return args[0].pow(args[1].coeff(0,0)); - - case MP2_SAFEPOWER: // SAFE_POWER, equivalent to (ifelse(x==0, 0, x^y)) - + + case MP2_SAFEPOWER: // SAFE_POWER, equivalent to (ifelse(x==0, 0, x^y)) + if (n != 2) { + SetError(err_code, "safe_power requires exactly two arguments", row, table_x[row] + 1, args.all_rows(), args.all_cols(), args.all_type_ints()); + return m; + } args = args.recycle_for_bin_op(); err_code = args.get_error_code(); switch (err_code) { @@ -1310,11 +1312,11 @@ class ExprEvaluator rows = m1.rows(); cols = m1.cols(); m = matrix::Zero(rows, cols); - for (int i; i < rows; i++) { - for (int j; j < cols; j++) { + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { m.coeffRef(i, j) = CppAD::CondExpEq( m1.coeff(i, j), - m2.coeff(i, j), + Type(0), Type(0), pow(m1.coeff(i, j), m2.coeff(i, j)) ); diff --git a/tests/testthat/test-binop.R b/tests/testthat/test-binop.R index abb618e4b..44ee112da 100644 --- a/tests/testthat/test-binop.R +++ b/tests/testthat/test-binop.R @@ -1,3 +1,5 @@ +library(macpan2); library(testthat); library(dplyr); library(tidyr); library(ggplot2) + test_that("elementwise binary operator executable specs match spec doc", { ## https://canmod.net/misc/elementwise_binary_operators times = BinaryOperator(`*`) @@ -83,3 +85,30 @@ test_that("elementwise binary operator executable specs match spec doc", { test_that("equivalent unary and binary minus operators give the same answers", { expect_equal(engine_eval(~-4), engine_eval(~0-4)) }) + +test_that("safe_power meets the requirements of #200", { + expect_equal( + matrix(0.1 ^ -5), + engine_eval(~safe_power(0.1, -5)) + ) + expect_equal( + matrix(0 ^ 5), + engine_eval(~safe_power(0, 5)) + ) + expect_equal( + matrix(c(1, 1, 0, 1, 1)), # != matrix((-2:2)^0), + engine_eval(~safe_power(-2:2, 0)) + ) + expect_equal( + matrix(rep(0, 5)), # != matrix(0^(-2:2)) + engine_eval(~safe_power(0, -2:2)) + ) + expect_error( + engine_eval(~safe_power(1)), + regexp = "safe_power requires exactly two arguments" + ) + expect_error( + engine_eval(~safe_power(1:3, 2:3)), + regexp = "The two operands do not have the same number of rows" + ) +})