Behavior of pow(x, y) when y is a variable with integer value

[perrydv]

Per Brad's email of two days ago, I'm starting a GitHub Discussion instead of using the mailing list to ask my question.

I get that handling of pow(x, y) is tricky because the calculations are not well defined if x < 0 and y is not an integer. CppAD provides pow(x, y) valid for x > 0 and pow_int(x, y) valid for any x when y is a constant integer.

I have use cases where CppAD taping is done by generic functions that call a lot of other functions that are themselves code-generated from end-user specifications. These are then used in many different ways.

What I would ideally like to have happen for pow(x, y), if both x and y are variables and x < 0, is that derivatives with respect to x are valid if y is integer-valued, NaN otherwise, and derivatives with respect to y are NaN. That would allow a user to get valid derivatives if only those with respect to x are sought and y is integer-valued, and to (correctly) get a NaN if invalid derivatives are sought. The current behavior is to abort.

If I follow, the implementation in pow_op.hpp is defined as having three outputs, the first two of which are internal (log(x) and y * log(x) in 0th order) for potential use in higher-order sweeps. If I make y a dynamic parameter, at least I can change its value via new_dynamic. If y has integer value and x <= 0, there are still NaNs from log(x) in the two internal outputs of pow_op, but they don't seem to propagate to cause a problem, as far as my small tests went. But in my generic use case, what I really need is for y to be a variable without triggering aborts based on particular values during sweeps.

Is there a good way to do this? Thanks for any ideas.

[perrydv]

A quick follow-up: I see that if I set check_for_nan to false, the abort behavior is avoided. But should I worry that the NaNs in the internal results of pow_op could still cause problems? It looks like first order reverse is correct but first order forward is NaN when it could be a valid result.

[bradbell]

You could write your own atomic function that computes pow(x, y); see
https://coin-or.github.io/CppAD/doc/atomic_three.htm
https://coin-or.github.io/CppAD/doc/atomic_three_example.htm

In fact, this is a good way to experiment with the design of a better computation.

[bradbell]

It might easier to do what you want using conditional expressions and a discrete function that returns the integer part of an AD value; see
https://coin-or.github.io/CppAD/doc/discrete.htm
https://coin-or.github.io/CppAD/doc/condexp.htm

You would compute the result two ways:

1. using logs
2. using pow_int and the discrete function.

Then, using the discrete function and a conditonal expression, you could choose pow_int when the AD value was equal to its integer part, otherwise use the log method.

[perrydv]

Thanks Brad. Those options make sense. It looks like the y of pow_int is baked in when the tape is recorded, so would that make it not possible to change later (I am looking at pow in utility/pow_int.hpp)? Would one approach be more efficient than the other? At one point in some efficiency experiments a while ago, it looked to me like there is substantial cost to the conditional expressions. Is that generally true? It could have been a quirk of whatever I was doing at the time.

[bradbell]

You are correct, one could only have a finite number of integer cases that would with conditional expressions. Perhaps you can make an atomic function that works for you. I would start with just function values, then add first order forward mode.

Otherwise, I do not know how soon I will be able to look into extending pow as you suggest above.

[bradbell]

My testing indicates that pow(x, y) for x < 0 and y integer works; see PowTestEight in
https://github.com/coin-or/CppAD/blob/master/test_more/general/pow.cpp

It is only the case where x = 0 and y is a positive integer that is a problem ?

[bradbell]

I have added an example that may solve your problem. Please take a look at
https://coin-or.github.io/CppAD/doc/pow_nan.cpp.htm

[perrydv]

Thanks for all of this! I'm thinking of the case where x and y are both variables. I modified your PowTestEight in https://github.com/perrydv/CppADtests/blob/main/PowTestEight_modified.cpp
to show the case. Execution aborts during first order forward. Based on your feedback yesterday, I am thinking that writing an atomic makes sense.

[bradbell]

I see the problem. In the representation z = exp( log(x) * y ), all the zero order coefficients get calculated correctly which sets the zero order coefficient for x to nan. Then when calcuating x^0 * y^1, which is used to compute the first order coefficient z^1, there is a multiply of nan * 0. It this case, it would make set ot use an absolute zero multiply (because the zero means we are not differentiating w.r.t. y), but in general, I do not want to use azmul for every multiply during forward mode.

[bradbell]

As an aside, you may find it useful, for small test cases, to use bin/trace.sh to turn on and off the tracing of the CppAD operation sequence (executing the script with no arguments gives a brief usage message). It will show exactly what operation a nan occurs at. For your test

    bin/trace.sh forward1.hpp 1

would should the operation. You need to use the include directory cppad.git/include for this to work (where cpapd.git is the directory where you cloned the repository).

[bradbell]

I recently did some improvement to the pow function; see discussion #93