Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Feature/issue 123 complex numbers with vars #789

Closed
wants to merge 141 commits into from
Closed

Feature/issue 123 complex numbers with vars #789

wants to merge 141 commits into from

Conversation

ChrisChiasson
Copy link
Contributor

@ChrisChiasson ChrisChiasson commented Mar 7, 2018 •
edited
Loading

Discourse use cases and spec discussion: http://discourse.mc-stan.org/t/complex-number-support-in-math-library-autodiff/3511

Submission Checklist

  • Run unit tests: ./runTests.py test/unit
  • Run cpplint: make cpplint
  • Declare copyright holder and open-source license: see below

Summary:

Initial attempt at adding complex number support to var, as referenced in issue #123, but extended further to do what I have laid out in the "Intended Effect" immediately following this.

Intended Effect:

Complex numbers work. For me, this means that it satisfies what I needed to do with complex numbers, which I described (after being prompted below) as:

Conceptually, and ignoring details, it is roughly equivalent to this: I have some tables of complex numbers from which I am indirectly generating linear combinations of some of the columns. Later, I am interpolating (via circular lerp) the resulting combination columns and integrating a function of the amplitudes. I have some checks for my integrations, and I will use multidimensional newton's method (requiring gradient and hessian) to minimize the error norm between the checks and the integrations as a function of the six coefficients of the aforementioned linear combinations and see if the result makes physical sense. My intention is to implement enough to accomplish this.

It also means that the PR passes the tests I have added for forward, reverse, and the special mixed mode test that uses newton's method and involves checking (the definition of) eigendecomposition of a rotation matrix while converging to the point where the merit function is minimized because the real and imaginary parts of the first returned eigenvalue are equal.

Further clarification for reviewers Friday 2018-07-27:
The paragraphs above mean that complex<var> and complex<fvar<T>> can be instantiated, added, subtracted, multiplied, divided, automatically differentiated, and integrated. As a side effect of getting the above working, the implementation was necessarily already very general and it was apparent that other functions would work. @bgoodri decided to take it further. He wrote additional tests (and I helped fix the resulting errors for) the list of C++ functions of complex numbers.

With respect to traits: without the additional definition for Eigen's scalar binary op trait, it would not be possible to have matrix operations on complex AD types, so that is how I came to the conclusion that it was necessary.

With respect to the fvar ctor that previously took a dummy parameter for SFINAE purposes: I ended up changing this to eliminate the dummy parameter due to compilation errors I was getting on a development machine, as previously discussed.

With respect to the previously existing var ctors that construct a var from the real part of a complex<double> and throw an error if the supplied complex<double>'s imaginary part is non-zero: I removed these because their original stated purpose was to check that Eigen didn't instantiate vars from complex numbers of this type. However, Eigen does not actually ever use those ctors, so this PR removes that dead code.

With respect to the large number of files touched: some of this is due to Jenkins overzealous autoformat.

How to Verify:

Run unit test

Side Effects:

One I can think of is that the interface of the complex class in stan math has a more flexible interface than the normal std::complex. This is because the binary operators use templates and ADL to do things like multiple a std::complex by a number without first casting the number to a var.

Documentation:

I added some inline documentation using doxygen syntax.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company): Oil States

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

bob-carpenter
bob-carpenter reviewed Mar 7, 2018
View reviewed changes
Copy link
Contributor

@bob-carpenter bob-carpenter left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for contributing. It'd be great to get some basic support for complex numbers. Is this intended to be something that eventually makes its way to the Stan language?

The first step needs to be a functional spec for what is being proposed---this doesn't look like it's still on the same course as the issue #123. Then a wider group can look at the spec and make sure it's going to work with al the various moving pieces of Stan.

We have very strict dependency requirements that this is violating. We do not allow files in the var path to include fvar, for example. I'm not sure where all of this is documented. I don't see anything in the math repo readme or Wiki.

I have a hard time imagining we're going to be able to get away with imaginary and real components being of different types. The problem is what happens when they combine with each other. But I take the point in the comments about type inference and templates---we run into that all the time. In general, you dont need to doc the language---assume the code reader is very familiar with C++, so that you only need to comment when you violate standard idioms.

@seantalts seantalts self-requested a review March 7, 2018 21:07
@seantalts
Copy link
Member

seantalts commented Mar 7, 2018 •
edited
Loading

Here is one page with some information sort of like this: https://github.com/stan-dev/stan/wiki/Code-Quality

Anyone know if we have a single page we can link new devs to? @syclik maybe?

Edit: Here are some pages that are decent to skim:
https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan
https://github.com/stan-dev/stan/wiki/Developer-process-overview
https://github.com/stan-dev/stan/wiki/Dev:-Git-Process
https://github.com/stan-dev/stan/wiki/Coding-Style-and-Idioms
https://github.com/stan-dev/stan/wiki/Code-Quality
https://github.com/stan-dev/stan/wiki/How-to-Write-Doxygen-Doc-Comments

@bob-carpenter
Copy link
Contributor

I see, what I was looking for (description of code layout and dependencies) is here:

https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

We need to prioritize getting that single getting-started page for new contributors.

@syclik
Copy link
Member

syclik commented Mar 7, 2018 via email

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Mar 7, 2018 via email

@seantalts
Copy link
Member

Don't worry! Things are a little complicated and we're still apparently working on our intro material, haha.

We've tried to keep var and fvar code completely separate so far. We have a mix directory that is for higher order code that would likely be using fvar<var> and fvar<fvar<var>>, but you only put code in there if you need to specialize for that case.

I think your code could go in prim, though we'd then take out the fvar and var forward declarations and the enable_if and lose a little bit of type safety there.

PS Do you have experience using forward declarations to reduce compile times? It's something I've been thinking about lately but haven't done before. Compile time is a big concern of ours.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Mar 7, 2018 via email

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Mar 8, 2018 via email

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Mar 8, 2018 •
edited
Loading

Our directories are organized as:

  • scal: int, double, var, fvar<T>
  • arr: std::vector, where Tis any Stan type
  • mat: Eigen::Matrix<T, R, C> where T is a scalar type and where R and C are the row and column int specifies determining if we have a vector (-1, 1), row vector (1, -1), or matrix (-1, -1).

That follows the data types in Stan of real and int for scalars, T[] for arrays, and vector, row_vector, and matrix for matrix types.

@bob-carpenter
Copy link
Contributor

Thansk for the clarification---you're absolutely right. Just a simple forward declaration isn't going to cause a bunch of code to be brought in as I was suggesting---it's not an include.

We still don't want to bring in the symbol because nothing will be able to bring in the declaration or definition.

@syclik
Copy link
Member

syclik commented Mar 8, 2018 via email

@syclik
Copy link
Member

syclik commented Mar 8, 2018

The way we have it, any templated code that can accept either primitives, var, or fvar<T> lives in prim. Template specializations then live in their respective folders.

I'll need to dig a little to figure out where it should live within scal.

@syclik
Copy link
Member

syclik commented Mar 8, 2018

I'm not sure. (I've also been sick with a headache, so not completely on... I'll take a look again tomorrow.)

@bob-carpenter, thoughts? I'd trust his judgement if he gets to it before I do.

@syclik
Copy link
Member

syclik commented Mar 8, 2018

I just looked back at the code a little more carefully. Take all I've said with a grain of salt. I may not have gotten all the nuance there. (I'm still trying to get around the fact that stan::math::internal::complex<T> is a subclass of std::complex<T>, but then we template specialize std::complex<stan::math::var> as a subclass of stan::math::internal::complex.)

So... I'm not sure exactly where it belongs. If we never need to template specialize for double, I'm almost tempted to suggest just specializing std::complex<var> and std::complex<fvar<T>> in their respective folders and not have anything live in prim. That would duplicate code -- based on the way you've handled the operators -- but I think it may be clearer? I'm not sure yet.

@syclik
Copy link
Member

syclik commented Mar 8, 2018

Yup. Makes sense.

@betanalpha
Copy link
Contributor

betanalpha commented Mar 8, 2018 via email

@betanalpha
Copy link
Contributor

betanalpha commented Mar 8, 2018 via email

@betanalpha
Copy link
Contributor

betanalpha commented Mar 8, 2018 via email

@seantalts
Copy link
Member

Since we need also need a spec before we can really evaluate a PR (as Bob mentioned), I went ahead and made a discourse thread asking for user input on applications. We can also use it for spec discussion (or we could make a new github issue for that, that could be better because we can summarize the comments at the top). Here's the link to the discourse thread: http://discourse.mc-stan.org/t/complex-number-support-in-math-library-autodiff/3511

@syclik
Copy link
Member

syclik commented Mar 8, 2018 via email

@seantalts
Copy link
Member

seantalts commented Mar 8, 2018 via email

@syclik
Copy link
Member

syclik commented Mar 8, 2018 via email

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Mar 9, 2018 via email

@seantalts
Copy link
Member

I think we're asking you to document what you're aiming to do before we merge your PR, haha.

@seantalts
Copy link
Member

We don't have ay specific format, just something that lets us know what operations you want to work, what's not supported, that kind of thing. If you're dead set on coding before getting feedback on a spec, your tests might function as a sort of very detailed and functional spec, which I think would be great.

@syclik
Copy link
Member

syclik commented Mar 27, 2018 via email

@seantalts
Copy link
Member

For Jenkins, you are getting cpplint errors - check out the make target for cpplint and see if you can run a similar command on your modified files, or just use the Jenkins output at the link (it should tell you which files are problematic and how to fix the issues).

For Travis, I'm not sure where you're seeing what you're seeing, but the errors I see immediately at the Travis job link are that you're using C++17 stuff?

./stan/math/cplx/complex.hpp:29:1: warning: inline variables are a C++17 extension [-Wc++17-extensions]
inline constexpr bool is_cplx_v = is_cplx<T>::value;

@seantalts
Copy link
Member

I misunderstood your question! Sorry about that.

I have had bunches of trouble with Ubuntu's C compiler packages when trying to use C++11. I have no idea why - seems fairly basic, but I am no Ubuntu wizard.

If you're trying to duplicate the travis environment, (which is Ubuntu 14.04, you can start with 14.04 and see how the .travis.yml file specified additional PPLs and packages to install: https://github.com/stan-dev/math/blob/develop/.travis.yml

So I think if you want clang, you want these additional PPLs: https://github.com/stan-dev/math/blob/develop/.travis.yml#L22

And then to apt-get install these packages: https://github.com/stan-dev/math/blob/develop/.travis.yml#L23

@seantalts
Copy link
Member

Also just want to confirm that you should indeed be using make/local for these flags (though in our travis config we just set environment variables). You can see how Jenkins does this in the Jenkinsfile

@syclik
Copy link
Member

syclik commented Jul 26, 2018

A few comments:

  1. The EXPECT_* macros have the first value as the expected value and the second as the actual value. So for line 39, instead of EXPECT_NEAR(f_x[pos++], 1, TOL);, it should be EXPECT_NEAR(1, f_x[pos++], TOL);
  2. For comparison against numbers that aren't supposed to be close to 0, use EXPECT_FLOAT_EQ() (or better yet, EXPECT_DOUBLE_EQ()... Google Test has come a long way since we've started).
  3. As @bgoodri said, relax TOL for comparison against 0.

@syclik
Copy link
Member

syclik commented Jul 26, 2018

Some relevant information:

@syclik
Copy link
Member

syclik commented Jul 27, 2018

@ChrisChiasson, as far as I can tell, this isn't ready for review. I still don't know what you're trying to do!

Can you write up a spec for what you're trying to accomplish and put it in the associated issue? If you have the ability, updating the first comment would be great.

Then, can you write some notes to guide the reviewer through this implementation against that spec that you describe? It'd help to know some of the major design decisions and whether you've implemented the whole spec or part of it and there's more remaining.

@syclik
Copy link
Member

syclik commented Jul 27, 2018 via email

@bgoodri
Copy link
Contributor

bgoodri commented Jul 27, 2018 •
edited
Loading

The use of the PR should be clear. If Stan Math <= 2.18 is an autodiff library for real numbers, then with this PR, Stan Math becomes an autodiff library for complex numbers conceptualized as a pair of var (or fvar<T>) one for the real part and one for the imaginary part. If z = a + b * i is a std::complex<stan::math::var>, then stan::math::foo(z) should create valid nodes on the autodiff tree, where foo is any function defined by the C++14 standard, such as log.

There are several situations in which you would want to autodiff functions involving complex numbers, but none of them are actually implemented in this PR, which is fine because this PR is just the framework for them. Future PRs could implement FFTs, various Eigen functions that utilize complex numbers internally, PDFs for complex random variables, etc. We should not wait on this PR for things like that.

@syclik your questions are largely about the implementation decisions. I don't fully understand why myself, but I do know that everything before this PR that pertained to complex numbers either would not compile or would segfault as soon as a complex var is instantiated. The main challenge is that the default constructor for a std::complex<T> is only well-defined when T is float, double, long double, etc. and it segfaults when T is a var. So, the guiding principle of the stuff in the cplx subfolder is to circumvent calling the default constructor, and it is quite possible that this PR is the only (known) feasible way to generalize Stan Math to be an autodiff library for complex numbers.

The second principle is that functions in Stan Math under rev and fwd need to get called on the real and imaginary parts of a std::complex<stan::math::var>. So, there is a bunch of stuff to ensure that happens, as opposed than calling functions in std or some other namespace.

But there are not many functions whose signature is a std::complex<stan::math::var> or similar. There are a few because even though C++11 defines a bunch of functions involving complex numbers that the standard library must implement, the implementations of them were not that good in earlier versions of stdlibc++ and libc++. With g++-4.9, there are several functions in the standard library that call other functions in the standard library but prefix them with std:: instead of allowing specializations to take precedence. So, there is some confusing stuff where we avoid calling functions in the standard libraries and basically reimplement what the standard libraries do or should do for those functions.

@syclik
Copy link
Member

syclik commented Jul 28, 2018

@syclik first post updated per your comments

Awesome. Thanks!

The use of the PR should be clear. If Stan Math <= 2.18 is an autodiff library for real numbers, then with this PR, Stan Math becomes an autodiff library for complex numbers conceptualized as a pair of var (or fvar) one for the real part and one for the imaginary part. If z = a + b * i is a std::complexstan::math::var, then stan::math::foo(z) should create valid nodes on the autodiff tree, where foo is any function defined by the C++14 standard, such as log.

There are several situations in which you would want to autodiff functions involving complex numbers, but none of them are actually implemented in this PR, which is fine because this PR is just the framework for them. Future PRs could implement FFTs, various Eigen functions that utilize complex numbers internally, PDFs for complex random variables, etc. We should not wait on this PR for things like that.

Thanks! That's the clearest I've seen it written out. Just to be clear, we never want to take a derivative of a std::complex<stan::math::var>, correct? We only want to take derivatives of the real or imaginary part (but never both)?

Thanks both for clarifying the issue. That explains what we want.

Time to get to the implementation and whether this is a decent implementation that gets us there. I'll start reviewing the code.

@syclik
Copy link
Member

syclik commented Aug 2, 2018

Mind merging develop back in? It looks like some of the diff isn't related to the issue at hand and it's a little distracting.

Since i is just a constant, the derivative of a complex number (which is actually a stupid name because it isn't plural) is the same as the derivative of any other sum of two terms: x' + constant * y', where constant is the value i. To me this means the derivative of both components is taken.

We're not using the same terminology. If x = f(theta) and y = g(theta), then when you say you want the derivative for both components, what I understand that to mean is to have both:

  1. dx / dtheta = f'(theta)
  2. dy / dtheta = g'(theta)

Instead, what you're returning is just one term (when theta is scalar):

  1. d (x + constant * y) / dtheta = f'(theta) + constant * g'(theta)

I just don't read that to mean the "derivative of both components is taken"; hope that clears up a lot of the confusion in this thread and on discourse. I now understand that what you mean and how that plays out, so it's a lot clearer.

I'll try to review this, but I may have a lot of questions for you.

syclik
syclik reviewed Aug 2, 2018
View reviewed changes
stan/math/fwd/scal/fun/is_inf.hpp
@@ -20,6 +20,12 @@ inline int is_inf(const fvar<T>& x) {
return is_inf(x.val());
}

// forwarding for ADL
template <class T>
inline auto isinf(const fvar<T>& a) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A bunch of questions:

  1. why is this function here? (this is isinf and not is_inf)
  2. why is it declared and defined multiple times? (in other files too)
  3. why is isinf in the stan::math namespace? I'm guessing you need it in the std namespace?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see how we've handled it elsewhere. Is there a reason we can't do the same here? By that, I mean specialize isnan and isinf for stan::math::fvar<T> and include it in its own header file called stan/math/fwd/core/std_isnan.hpp and std_isinf?
See:
stan/math/fwd/core/std_numeric_limits.hpp
stan/math/rev/core/std_isnan.hpp
stan/math/rev/core/std_numeric_limits.hpp
stan/math/rev/core/std_isinf.hpp

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you clarify? What other packages? I didn't see anywhere in this pull request that this is meant to work with other packages.

(I just found a stackoverflow post on isnan outside of any namespace; is this something you've run into or is this hypothetical?)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we're way past that; std::complex for our types are unspecified. We should do what we can verify is correct, is maintainable, and makes sense.

@syclik
Copy link
Member

syclik commented Aug 2, 2018

@ChrisChiasson, did you go through all the implementation notes here: https://en.cppreference.com/w/cpp/numeric/complex? (that's a question... I haven't verified that you did or didn't)

There's no discussion in this pull request regarding what parts of it you have or haven't implemented.

@seantalts seantalts removed their request for review August 2, 2018 14:56
@syclik
Copy link
Member

syclik commented Aug 2, 2018

@ChrisChiasson, maybe this would explain it?

#include <stan/math.hpp>
#include <iostream>


int main() {
  std::cout << "compute dx[0] / dtheta" << std::endl;

  std::vector<stan::math::var> x;
  std::vector<double> gradient;

  stan::math::var theta = 0.5;
  std::vector<stan::math::var> vars;
  vars.push_back(theta);

  x.push_back(sin(theta));
  x.push_back(cos(theta));

  std::cout << x[0] << std::endl
            << x[1] << std::endl
            << std::endl;

  x[0].grad(vars, gradient);

  std::cout << "gradient.size() " << gradient.size() << std::endl
            << "dx[0] / dtheta = " << gradient[0] << std::endl;

  // clear everything
  stan::math::recover_memory();
  x.clear();
  vars.clear();

  // d x[1] / dtheta
  std::cout << "compute dx[1] / dtheta" << std::endl;

  theta = 0.5;
  vars.push_back(theta);
  x.push_back(sin(theta));
  x.push_back(cos(theta));

  std::cout << x[0] << std::endl
            << x[1] << std::endl
            << std::endl;

  x[1].grad(vars, gradient);

  std::cout << "gradient.size() " << gradient.size() << std::endl
            << "dx[1] / dtheta = " << gradient[0] << std::endl
            << std::endl;

  // clear everything
  stan::math::recover_memory();
  x.clear();
  vars.clear();

  // compute both d x[0] / dtheta and d x[1] / dtheta simultaneously
  std::cout << "can't compute both dx[0] / dtheta and dx[1] / dtheta simultaneously" << std::endl;

  return 0;
}

@syclik
Copy link
Member

syclik commented Aug 2, 2018

The entire set of commits was rebased onto develop right before my review request.

Ok. Can you revert the changes you've made to the makefiles?

@syclik
Copy link
Member

syclik commented Aug 2, 2018

Is there some part of it you would like me to comment on?

Sure. All of it? Did you satisfy all the requirements? How did you validate that std::complex<stan::math::___> is a LiteralType? Did you check that the the real and imaginary parts are sitting next to each other in memory?

There's a lot there... I don't know how much you walked through it or didn't.

@syclik
Copy link
Member

syclik commented Aug 2, 2018

re: make. ok. I have no idea what GitHub is doing.

@syclik
Copy link
Member

syclik commented Aug 2, 2018

@ChrisChiasson, yes, I read that.

Questions that would help everyone else understand what you're doing:

  • what parts of the (soft) requirements are you implementing?
  • what parts are you explicitly not implementing?
  • why? (to both)
  • did you add tests that checked those things?

In a pragmatic sense, things that are undefined are ok as long as we know that the compiler does the right thing (it looks like they start allowing it for C++20). If we have tests in place, at least we can validate that we're ok even if we're in undefined behavior land.

@syclik
Copy link
Member

syclik commented Aug 2, 2018

re: contiguous. Awesome. I didn't realize that was going to be implemented by the standard library. It'd just help line up these things.

re: limitation of grad. Hopefully that clears it up for you why I was asking for what you were attempting to do. That's not a limitation in our design: it's a limitation of reverse-mode automatic differentiation when it can't be retaped. Also, I don't understand what it means to compute .grad() on std::complex -- I'm not sure what that would return. Anyway, I understand that it's not what we're trying to do, so let's focus on helping me (or any other reviewer) understand what you're accomplishing and how.

syclik
syclik reviewed Aug 6, 2018
View reviewed changes
stan/math/fwd/scal/fun/is_inf.hpp
@@ -20,6 +20,12 @@ inline int is_inf(const fvar<T>& x) {
return is_inf(x.val());
}

// forwarding for ADL
template <class T>
inline auto isinf(const fvar<T>& a) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see how we've handled it elsewhere. Is there a reason we can't do the same here? By that, I mean specialize isnan and isinf for stan::math::fvar<T> and include it in its own header file called stan/math/fwd/core/std_isnan.hpp and std_isinf?
See:
stan/math/fwd/core/std_numeric_limits.hpp
stan/math/rev/core/std_isnan.hpp
stan/math/rev/core/std_numeric_limits.hpp
stan/math/rev/core/std_isinf.hpp

stan/math/fwd/scal/fun/proj.hpp
@@ -0,0 +1,28 @@
#ifndef STAN_MATH_FWD_SCAL_FUN_PROJ_HPP
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just following our (not very good) convention, can you rename this file std_proj.hpp?

stan/math/prim/cplx/complex.hpp
namespace math {

/**
* Complex class that forwards the interface of std::complex, brining
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is meant here? I still don't understand what "forwards" means in the way it's used here.

stan/math/prim/cplx/complex.hpp
* division operations that the base std::complex doesn't have
* defined usably for var, and, since stan::math::zeroing<var> zero
* initializes itself, will also correctly work with the remaining
* algorithms in std::complex<T> that require the expression T() to
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clarification? The algorithms expect T() to initialize a variable to 0?

stan/math/prim/scal/meta/rm_zeroing.hpp
* T. This is the removal version of the trait.
*/
template <class T>
struct rm_zeroing<math::zeroing<T>> {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove math:: in front of zeroing

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

at that line in the file, the math sub-namespace is not open, just the outer stan, so if I did that, it would cause a compile error

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I didn't catch that on the first read. Forgot that our meta is in the Stan namespace.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I didn't catch that on the first read. Forgot that our meta is in the Stan namespace.

stan/math/prim/cplx/complex.hpp
* be 0.
*/
template <class T>
struct complex : std::complex<zeroing<T>> {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the motivation for having stan::math::complex<T>? Shouldn't we just be using std::complex<T>?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we did that, we would create an inheritance loop. The way we force things to work is by making a template specialization for std::complex<var> inherit from stan::math::complex<var>, which in turn inherits from std::complex<zeroing<var>>. We wouldn't be able to do this if the starting and ending classes were the same.

@syclik
Copy link
Member

syclik commented Aug 6, 2018

@ChrisChiasson, I've been staring at the code. I just made some minor-ish comments (been sitting there as I tried to make heads / tails of the code).

Some bigger questions that don't fit into a line-by-line commentary:

  1. Why is there a stan::math::complex<T> type defined?
  2. Instead of having the zeroing construction, would it be possible to just write the instantiation so it calls the constructor with 0? Is it difficult to completely template instantiate the class? Are you trying to avoid that?
  3. Can you explain the side effect a little more?

Sorry it's taken a while. Here's what would help me (or anyone else) review it:

  • a guide to how to navigate your code.
  • any design decisions you made: why zeroing?
  • what functions are implemented in the std namespace and why
  • what functions are implemented in math? (just to have a list)

I haven't walked through the tests yet. It'd help to know what you were trying to test. If not, it'll take a while to walk through and understand what you're trying to test.

@syclik
Copy link
Member

syclik commented Aug 7, 2018

@ChrisChiasson, I probably should have asked this earlier on. It just dawned on me.

In your experimentation, did you ever write a template specialization for any of these?

  • std::complex<stan::math::var>
  • std::complex<stan::math::fvar<double>>
  • std::complex<stan::math::fvar<stan::math::var>>
  • std::complex<stan::math::fvar<stan::math::fvar<double>>>
  • std::complex<stan::math::fvar<stan::math::fvar<stan::math::var>>>

Would it be possible to write an explicit template specialization for just the constructor? Here's a short example along the lines of what I'm thinking (but I haven't done it for complex):

#include <stan/math.hpp>
#include <iostream>

template <typename T = double>
class foo {
 public:
  foo(T x = T(), T y = T())
      : x_(x), y_(y) {
  }
  std::string print() {
    std::stringstream m;
    m << "(" << x_ << ", " << y_ << ")" << std::endl;
    return m.str();
  }
  T x_, y_;
};

template<> foo<stan::math::var>::foo(stan::math::var x, stan::math::var y) {
  if (is_uninitialized(x))
    x_ = 0.0;
  if (is_uninitialized(y))
    y_ = 0.0;
}

int main() {
  foo<double> d;
  foo<stan::math::var> v;
  std::cout << "d = " << d.print() << std::endl;
  std::cout << "v = " << v.print() << std::endl;
  return 0;
}

If we could just write a constructor for each, that would eliminate the need for zeroing, the template metaprograms (e.g. rm_zeroing), and simplify a bunch of it? I'm not sure if you've already explored this and found that it wasn't possible.

@syclik
Copy link
Member

syclik commented Aug 7, 2018 via email

@syclik
Copy link
Member

syclik commented Aug 8, 2018

Just to take inventory:

  1. the default constructor (which we can template specialize): https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex#L129
  2. operator= for assigning a real to a complex (which we should be able to specialize?): https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex#L235
  3. operator== for complex, real. Isn't this ok? https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex#L463
  4. operator!= for the same reasons: https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex#L481
  5. __complex_abs. I don't think there's anything wrong with this implementation. https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex#L593

There are another half dozen places where it pops up. I don't see a problem with them, do you?

@syclik
Copy link
Member

syclik commented Aug 8, 2018

@ChrisChiasson, I searched through your tests. What I didn't find were basic tests for std::complex<stan::math::var> and all the other types. Can you either point me to them or can you create the tests for them?

This would be a requirement for accepting this sort of code. What I mean by tests:

  • test that you can instantiate std::complex<T> where T is {stan::math::var, stan::math::fvar<double>, stan::math::fvar<stan::math::var>}. You could do fvar<fvar<var>>, but I think that'd cover it.
  • test for each of the operators: (*, /, +, =, !=, *=, +=, etc.); it looks like we should be testing with complex arguments and real arguments (based on what you pointed out with the_Tp())
  • if you're so inclined, you should test the few functions that have the _Tp() in the implementation.

Why here and not in something like std::vector<T>? It's because std::vector<T> is specified for custom types and we obey all the constructs. Here, the implementation may change and the only thing we have to know that it's still good is if we have tests.

I have a feeling that once we create these tests, we can replace the implementation with something a lot simpler: getting rid of zeroing, rm_zeroing, etc. And then we can merge easily.

@syclik
Copy link
Member

syclik commented Aug 8, 2018 via email

@syclik syclik mentioned this pull request Aug 9, 2018
Closed
@ChrisChiasson ChrisChiasson deleted the feature/issue-123-complex-numbers-with-vars branch August 9, 2018 18:02
@syclik syclik mentioned this pull request Sep 17, 2018
Closed
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@syclik syclik syclik left review comments

@bob-carpenter bob-carpenter bob-carpenter left review comments

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
10 participants
@ChrisChiasson @seantalts @bob-carpenter @syclik @betanalpha @bgoodri @bbbales2 @stan-buildbot @rok-cesnovar @SteveBronder