Fully functional std::complex specializations and overloads

Well, I intended all those tests to be evaluated across a bunch of random matrices in order to find out whether there are places where the accuracy gets unacceptable. But people hated non-deterministic unit tests.

lol, sorry :) The reason I tested them in mine is that I didn’t know exactly what Eigen was doing, specifically if there were any numbers initialized through identities they exploited in code without calculation, possibly breaking the adjoint flow. Might not be possible anyway; didn’t know / wanted to program defensively.

Hrm, if a is a reference, maybe the result of decltype(a) could be post-processed via decay or remove_reference?

Are the tolerances treated as absolute residuals (abs(left-right)<tol), as relative ones (abs(left-right)<tol*max(abs(left),abs(right))), or combined like (abs(left-right)<tol*max(1.0,max(abs(left),abs(right))))? I usually use a combined type of check to handle cases with large and small magnitudes, with tol around half the precision of a double (about 1.0e-8, or more closely pow(2.,-53./2.)). If it still isn’t loose enough, I provide a multiplier on top of tol. A combined approach should allow the tol to be a fairly small number unless the finite difference is in a really bad spot.

Just think how far we’ve come from when people were debating about whether or not it was even possible to have stan do this at all…

I’m reading your code more closely today before swapping out my implementation for yours in my application at work.

Right now, the inheritance structure I’m reading is:
std::complex<stan::math::var> inherits from stan::math::complex_base<stan::math::var>. The operations of stan::math::complex_base statically downcast their return references back to std::complex<stan::math::var>.

For reference, my inheritance structure is:
std::complex<stan::math::var> inherits from stan::math::complex<stan::math::var>, which in turn inherits from std::complex<stan::math::zeroing<stan::math::var>>>. There is also an operator std::complex() defined on stan::math::complex, which implicitly statically downcasts to std::complex for both base classes. (The extra level of inheritance together with the implicit downcasting delegates all operations back to the standard library as if they were operating on a zeroing version of stan’s types. This is why there is so little code, excluding the Jenkins autoformat misconfiguration that had touched a lot of files.)

While trying to test your code:
I was trying to build the libraries. I ended up partially following a page I wrote in the math wiki a while back, along with an example Jenkins test run. make test-headers fails with the error below. ./runTests.py fails because due to a missing Google Test library. I’m on commit 2c746c2cc9d19feb036dae474444b86d940be210 (HEAD -> feature/0123-complex-var, origin/feature/0123-complex-var). Edit: Forgot to add that I’m on Windows with MSYS2 x86_64 using GCC 9.1.

g++ -std=c++1y -m64 -D_REENTRANT -Wall -Wno-unused-function -Wno-uninitialized -Wno-unused-but-set-variable -Wno-unused-variable -Wno-sign-compare -Wno-unused-local-typedefs      -I lib/tbb_2019_U8/include -O3  -I . -I lib/eigen_3.3.3 -I lib/boost_1.69.0 -I lib/sundials_4.1.0/include  -D_USE_MATH_DEFINES  -DBOOST_DISABLE_ASSERTS      -c -std=c++1y -m64 -D_REENTRANT -Wall -Wno-unused-function -Wno-uninitialized -Wno-unused-but-set-variable -Wno-unused-variable -Wno-sign-compare -Wno-unused-local-typedefs      -I lib/tbb_2019_U8/include -O3  -I . -I lib/eigen_3.3.3 -I lib/boost_1.69.0 -I lib/sundials_4.1.0/include -O0 -include stan/math/mix/mat/fun/typedefs.hpp test/dummy.cpp -o nul
In file included from ./stan/math/rev/core.hpp:49,
                 from ./stan/math/mix/mat/fun/typedefs.hpp:6,
                 from <command-line>:
./stan/math/rev/core/std_complex.hpp: In instantiation of 'std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&) [with T = stan::math::var]':
./stan/math/rev/core/std_complex.hpp:693:48:   required from 'std::complex<_Tp> stan::math::complex_asinh(const std::complex<_Tp>&) [with V = stan::math::var]'
./stan/math/rev/core/std_complex.hpp:1537:25:   required from here
./stan/math/rev/core/std_complex.hpp:229:23: error: no matching function for call to 'value_of_rec(stan::math::complex_base<stan::math::var>::value_type)'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |           ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note: candidate: 'template<class T> std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&)'
  228 | inline std::complex<double> value_of_rec(const std::complex<T>& z) {
      |                             ^~~~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note:   template argument deduction/substitution failed:
./stan/math/rev/core/std_complex.hpp:229:23: note:   'stan::math::complex_base<stan::math::var>::value_type' {aka 'stan::math::var'} is not derived from 'const std::complex<_Tp>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |           ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:229:47: error: no matching function for call to 'value_of_rec(stan::math::complex_base<stan::math::var>::value_type)'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                   ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note: candidate: 'template<class T> std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&)'
  228 | inline std::complex<double> value_of_rec(const std::complex<T>& z) {
      |                             ^~~~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note:   template argument deduction/substitution failed:
./stan/math/rev/core/std_complex.hpp:229:47: note:   'stan::math::complex_base<stan::math::var>::value_type' {aka 'stan::math::var'} is not derived from 'const std::complex<_Tp>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                   ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:229:57: error: could not convert '{<expression error>, <expression error>}' from '<brace-enclosed initializer list>' to 'std::complex<double>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                                         ^
      |                                                         |
      |                                                         <brace-enclosed initializer list>

Thanks. I hadn’t tested the headers—those are always easy to fix with more includes, so I’ll fix them and check in a cleaner version. Aside from failing the header tests, were there other problems?

The reason I was reluctant to delegate back to the standard library is that instantation of std::complex<T> for anything other than a primitive type is undefined behavior.

I just used a simple CRTP base class to avoid shared code between forward and reverse std::complex implementations, figuring that cutting-and-pasting a whole interface wouldn’t get by code review.
And I wanted to support mixed operations of things like std::complex<double> and std::complex<var> without promotion. Eventually, if we care about performance, we’ll need to reimplement reverse mode with analytic derivatives, not just autodiff.

I didn’t find the need to zero anything anywhere other than in the std::complex constructors.

I’ll give it a shot. I just took the easy way out because I was wrestling with so many type errors on the way to getting everyting to compile.

The errors are relative unless one of the values is very small at which point it kicks over to absolute errors. The tolerance is set around 1e-8 for gradients, but by the time you get to Hessians with finite differences (and they’re based on three finite diffs, not just one in a band, which is partly why they’re so slow), things are much less stable arithmetically. The tolerance there is 1e-3 and often has to go down to 5e-3 or 1e-2. Longer term, I think the right thing to do is instead do finite differences of gradients, since the gradients are being independently tested. Same for third order derivatives, but those are so often zero that it’s not nearly as much of a problem for testing as Hessians.

We could also use more of the algebraic tests if people are concerned about precision. I just figured if Hessians were right up to 1e-2 they were probably actually right. That’s obviously not a watertight line of reasoning.

I don’t like that they can fail randomly and I don’t like that the coverage is random. I’d rather have someone think out the edge cases and test those.

I’m one of those people. The reasons are that They’re unruly in that they can fail randomly, especially with finite diffs. Isn’t it also hard to actually get edge tests like highly ill-conditioned matrices randomly?

I don’t like that the test coverage is a measure zero subset of the parameter space it is going to get used for with HMC. In general, our unit tests primarily test whether stuff compiles and whether the analytical derivatives are logically correct, but there is very little testing of whether we are getting the reasonable numerical values considering that numerical accuracy is critical to HMC. Checking that you get reasonable numerical values for a small matrix often with integer values does not tell you whether you are going to get reasonable numerical values when plugged into an actual Stan program. We have some experience as to what sort of things tend to go badly in Stan programs, but no accumulated experience as to what is going to happen when you extend out into the complex plane.

To test “edge” cases, you usually have to construct things manually, but it is not as if you need an edge case like a singular matrix for things to break down numerically. Often it can be just due to having big values in some cells and small values in other cells, particularly if the matrix is large. To find things like that out, it is often good to use random matrices, which is what Eigen largely does in its unit tests with basic scalar types.

What about using randoms, but if it hits an error, making it dump the input to enable devs to add a new test with exactly those inputs that gets run every time? Over the years, it might accumulate some really good test cases just by chance.

Hrm, sounds reasonable to me.

+1

Hrm, although I would note in the first instance that the goal of both implementations is already in conflict with the ‘rule’: “Specializing the template std::complex for any type other than float, double, and long double is unspecified.”

Of course, I never follow the rules anyway. Like using auto all the time with Eigen… Eh, to each his or her own. I’m just happy someone in a position to do something about it is adding stuff I need to stan.

I got further more recently by checking out master from develop, building a core rev test to generate the libs as a side effect, and then switching back to the complex branch. To integrate with my app, I had to figure out MinGW with both static and shared libs (tbb)… (side note, rpath does not work on Windows, and TBB likes to spawn 3 threads just for a single threaded unit test… ::arches eyebrow::).

I’ll check when I get in to work tomorrow, but it was telling me that pow(var,int) (I think that’s what it said) was ambiguous, an error that didn’t show up when using my implementation. I don’t have clang++ fully set up at work because the link step fails in MinGW, but I did copy down and email to myself the list of candidate functions from it (g++ doesn’t dump that info):


../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2568:26: note: candidate function
inline std::complex<var> pow(const var& lhs, const std::complex<double>& rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2585:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, int rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/scal/fun/pow.hpp:124:12: note: candidate function
inline var pow(const var& base, double exponent) {
           ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2589:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, double rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/scal/fun/pow.hpp:108:12: note: candidate function
inline var pow(const var& base, const var& exponent) {
           ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2593:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, const var& rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2597:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs,
                         ^

I dont want to derail this topic but not sure what you mean by “not working”? In order to do dynamic linking on Windows TBB needs to be in the users or systems PATH variable.

For Stan Math we let the TBB threadpool handle the number of spawned threads with their default policy. You can override that by calling stan::math::init_threadpool_tbb();. If STAN_NUM_THREADS is set, that call will spawn the number of threads set in that variable, otherwise it will spawn 1 threads. Cmdstan calls that at the start of each model.

rpath, if I understood it correctly, is different from the PATH environment variable. I believe you’re correct about PATH (in that it is one of the ways to solve the problem). I had only mentioned rpath because it is in the instructions for how to compile a stan math program in the readme.md.

It was pow(var,int). The beginning of the error snippet was as follows:

src/Nautical/MotionTables/MotionTables.hpp:55:16: error: call of overloaded 'pow(stan::math::var, int)' is ambiguous
   55 |        ess+=pow(ERSAft-check,2);
      |             ~~~^~~~~~~~~~~~~~~~
compilation terminated due to -Wfatal-errors.
In file included from src/main.cpp:9:
./src/Nautical/MotionTables/MotionTables.hpp:55:13: fatal error: call to 'pow' is ambiguous
       ess+=pow(ERSAft-check,2);

I have been wondering, why are you calling it CRTP? The pattern I see is

template <>
class complex<stan::math::var>
    : public stan::math::complex_base<stan::math::var> {

I agree with that and consider it worthwhile. I would only note the approaches are not mutually exclusive. In my implementation, sometimes we selectively avoided delegating back to the standard library, especially with g++ (as you mentioned), by writing declarations with more specific types to give them higher priority in template selection and overload deduction.

A finite set of random tests is still measure zero. How about rather than having unit tests, we run a million trials of something random under development? That way, we don’t have random inexplicable failures of unit tests unrelated to the branch being developed.

I called it CRTP-like in the doc. It would be an instance of CRTP if the template type were std::complex rather than T; but it’s used just like the CRTP in terms of casting to the subclass.

I’ll check out what’s going on. I may not have been testing that case.

I can try to fire up your PR and see if it passes all of my tests. I thought this was going to be a lot easier than it is.

Also, as of now, everything’s broken while I try to sort out the ambiguities in the header tests.

Indeed, pow(var, int) is also ambiguous for me with clang++:

  auto fx = pow(x, 2);
            ^~~
./stan/math/rev/core/std_complex.hpp:2599:26: note: candidate function
inline std::complex<var> pow(const var& lhs, const std::complex<double>& rhs) {
                         ^
./stan/math/rev/core/std_complex.hpp:2616:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, int rhs) {
                         ^
./stan/math/rev/scal/fun/pow.hpp:124:12: note: candidate function
inline var pow(const var& base, double exponent) {
           ^
./stan/math/rev/core/std_complex.hpp:2620:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, double rhs) {
                         ^
./stan/math/rev/scal/fun/pow.hpp:108:12: note: candidate function
inline var pow(const var& base, const var& exponent) {
           ^
./stan/math/rev/core/std_complex.hpp:2624:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, const var& rhs) {
                         ^
./stan/math/rev/core/std_complex.hpp:2628:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs,

Looks like I have some template wrangling to do to fix the calls. This one can be fixed by defining pow(var, int).

I’ve pushed fixes so that the header tests now pass and pow(int, var) and pow(var, int) work.

The unit tests are still failing to compile in various places due to ambiguities.

It would be interesting to see what would happen testing my PR. In many ways, I envy the testing framework you guys have now, because of the finite difference checks and the different type instantiations. Back when I did my implementation, I was mortally afraid of results being wrong (especially given what I’m using it for), so I wrote as little code as possible to leave less surface area for errors.
TL;DR: In my day, we walked to school in the snow uphill both ways. We ate dirt, and we were thankful.

Understood
I’ll be back in the office Monday & will try to check. Also, I may be slow to respond next week due to an industry meeting I have to attend. Apologies in advance.

1 Like

For me, because of the choice to delegate back to the standard library implementation, and because of g++ using _Tp() everywhere to represent zero, there were many functions where I needed zeroing: real assignment (239), most of the comparator ops (465, 470, 483, 488), abs (589), sqrt (894), pow(1021). Clang wasn’t nearly as bad, at least in that one respect.

No worries. I’m also traveling and haven’t been doing a lot around this. I’ll be back in the office next week myself.

I updated to commit 0769d931462b861ae08cff6129d88b0d9b3556fa (HEAD -> feature/0123-complex-var, origin/feature/0123-complex-var). I get

./src/Nautical/RAO/RAO.hpp:80:78: error: no match for 'operator<<' (operand types are 'Eigen::Matrix<std::complex<stan::math::var>, 3, 3, 0, 3, 3>' and 'double')
   80 |   return (Eigen::Matrix<decltype(std::complex(f_*phi_x*theta_y*psi_z)),3,3>()<<
      |                  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~
   81 |          0., -a_*psi_z, b_*theta_y,
      |          ~~

Update: ran it through clang++ instead of g++ above. The error is the same as g++, but clang dumped more information on the failed overloads. It looks like it cannot find a way to construct a complex<var> from a double. (Edit note: I had said this backwards earlier today)

../eigen-git-mirror\Eigen/src/Core/DenseBase.h:312:31: note: candidate function not viable: no known conversion from 'double' to 'const Eigen::DenseBase<Eigen::Matrix<std::complex<stan::math::var>, 3, 3, 0, 3, 3> >::Scalar' (aka 'const std::complex<stan::math::var>') for 1st argument
    CommaInitializer<Derived> operator<< (const Scalar& s);

1 Like

I haven’t tried to overload for Matrix<T, 3, 3> as we only use dynamic matrices and vectors in Stan, not fixed size ones. I’ll see if I can get this to compile, but it’s not generally in scope for Stan. Very little of our matrix code is general enough to handle Matrix<T, 3, 3>.