add 'safe_power' engine function that sets 0^0 to 0

[bbolker]

it would be convenient to have an engine function that computes x^y but always returns 0 if x==0. (We have a special case where we are computing 0^0, where C++ quite reasonably returns NaN; we know we want the answer to be 0 in this case.) I took a crack at this on a branch (with CppAD::CondExpLt), but am so far failing (perhaps because of some issue with vectorization? Not sure ...)

If this is not something @stevencarlislewalker can figure out quickly I might ask on the TMB list ...

(attn: @jfree-man)

[stevencarlislewalker]

Thanks Ben. I'll checkout your branch.

[bbolker]

PS There seems to be a similar construction at line 130 of https://kaskr.github.io/adcomp/convenience_8hpp_source.html

[stevencarlislewalker]

Thanks! I'm following this here: https://www.coin-or.org/CppAD/Doc/condexp.htm

[bbolker]

New branch at https://github.com/canmod/macpan2/tree/safe_power2, PR

[stevencarlislewalker]

Thanks Ben. I think I did it a while ago, then needed to take little people to choir and then there are tests failing from other commits in the branch. Here's the pull request I was using: https://github.com/canmod/macpan2/pull/201 . But I haven't pushed my changes yet. In any case, this is what it will look like I think.

case MP2_SAFE_POWER: // SAFE_POWER
                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<Type>::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;

[stevencarlislewalker]

Looks like CondExpLe is more appropriate.

[bbolker]

Why? In this case I don't think we want to be more permissive (i.e. safe_power(x,y) should be NaN for x<0)? #202 is the correct PR, I've added your code to it ...

[stevencarlislewalker]

I now have it working with tests, and was about to check it in when I realized that R itself seems to be doing ifelse(y == 0, 1, x^y). I'll check in my changes to the branch as they are (i.e. ifelse(x == 0, 0, x^y)`) but let me know if you think we should switch.

[stevencarlislewalker]

More confusion. I've now asked the engine to do the case in the description of this ticket, and I get the following.

> engine_eval(~0^0)
     [,1]
[1,]    1

This is not a NaN. Anyways, I'm still going to check in to the branch, but we should continue this discussion.

[stevencarlislewalker]

Another example of R doing something that my current safe_power doesn't emulate.

> matrix(0^(-2:2))
     [,1]
[1,]  Inf
[2,]  Inf
[3,]    1
[4,]    0
[5,]    0

[stevencarlislewalker]

Here is the first version of the test of safe_power.