coin-or / CppAD

A C++ Algorithmic Differentiation Package: Home Page
https://cppad.readthedocs.io
Other
446 stars 94 forks source link

How to add self defined functions to AD<Base>? #147

Closed xuruilong100 closed 2 years ago

xuruilong100 commented 2 years ago

I am trying to add Gamma function in this way:

template <class Base>
inline AD<Base> gamma(const AD<Base> &x)
{   
    return x.gamma_me();
}

How to implement it? Which parts of source codes should I modify?

bradbell commented 2 years ago

Does the following work for you ?

template <class Float>
Float gamma(const Float& x)
{
    // your C++ code for the gamma function
    ...
   return result;
}
xuruilong100 commented 2 years ago

In fact, there are a lot of 'if-else' branches in the numerical algorithm of Gamma function, see boost's implement

https://www.boost.org/doc/libs/1_78_0/boost/math/special_functions/gamma.hpp

I am worried about that the code may be non-differentiable.

bradbell commented 2 years ago

The operation sequence will depend on the results of all the if branches in your program; see https://coin-or.github.io/CppAD/doc/compare_change.htm

It is more complicated to create a gamma function where you never have to retape. For an example see the atan2 example https://coin-or.github.io/CppAD/doc/condexp.htm#Atan2

It is even more complicated, but more efficient, to create an atomic gamma function; see https://coin-or.github.io/CppAD/doc/atomic_four_get_started.cpp.htm Note that to use this in derivative calculations, you will have to define the derivatives of the gamma function.

xuruilong100 commented 2 years ago

Great, thank you very much. I have to say I need more tutorials about atomic_four function.

bradbell commented 2 years ago

The important thing is to only implement the virtual functions that you need. For more examples see; https://coin-or.github.io/CppAD/doc/atomic_four_example.htm

xuruilong100 commented 2 years ago

Really helpful! And can you tell me what do the args mean, for instance select_y and call_id.

bradbell commented 2 years ago

See the documentation for the particular function below https://coin-or.github.io/CppAD/doc/atomic_four_define.htm

bradbell commented 2 years ago

If you are satisfied with this answer, would you please close this issue.