Bugfix/issue 1250 lgamma

Summary
The std::lgamma is not reentrant safe and has been shown to avoid performance scalability in threaded stan programs using map_rect-type  threading on macOS. This PR replaces this with a call to the reentrant safe variants with one of these two:

lgamma_r the platform is not __MINGW32__
boost::math::lgamma on MINGW32

Both versions are reentrant safe. The lgamma_r is a C library call which is only available upon defining -D_REENTRANT, but this is not working on Windows with MinGW32. This is why the boost version is used on Windows. However, the boost lgamma implementation currently uses as default policy to propagate its input arguments from double to long double. This leads a serious performance regression as the boost lgamma is unacceptably slow in comparison to the std::lgamma (see below). Changing the default policy of boost math functions to avoid the propagation to long double leads to speed of the same order as std::lgamma, see below.
The proposed log gamma implementations solve the performance problems seen with threading under macOS.
Tests
The makefile target test-math-dependencies has been extended to check in stan/math for use of std::lgamma. In the commit 6521f3f the reviewer can verify that a
make test-math-dependencies
will rightfully detect the use of std::lgamma which is put in this commit for testing purposes into the source.
The boost math lgamma implementation does overflow for arguments greater than 1E25 whenever the double propagation to long double is disabled. This is a bug in boost 1.69 and will be fixed past boost 1.70.0. This PR includes the fix from the boost math sources (https://github.com/boostorg/math/commit/fd07b121a619e0429e189edb4ff50f5b4e9efe58). In order to check that the boost version is fixed a test is included in
test/unit/math/prim/scal/fun/lgamma_test.cpp
to check for large arguments to return a finite value. The bug appeared for arguments greater than 1E25 (around the sqrt of the largest double of 1E308). Other functions from boost math which use the stan::math::boost_policy_t (see below) are also tested to work for large arguments.
Side Effects
Change to the default boost_policy_t
We change the default boost_policy_t within stan math to avoid the double to long double propagation (this being considered to be the new default in future boost versions due to the performance regression caused, see https://github.com/boostorg/math/issues/241). As a side effect all functions from boost math which use this policy are affected. These are:
stan/math/prim/scal/fun/rising_factorial.hpp
stan/math/prim/scal/fun/falling_factorial.hpp
stan/math/prim/scal/fun/digamma.hpp
stan/math/prim/scal/fun/lgamma.hpp

To verify that things are still fine, I have run a negative binomial model which uses the digamma function and results are just the same with and without the overpropagation - with the difference of getting things ~2x faster on this branch.
Change to the C++ compiler flags: adding -D_REENTRANT
Now the  makefiles define -D_REENTRANT since otherwise cmath will not expose lgamma_r. The flag is set by default on all platforms in our makefiles (CXXFLAGS_OS).
Performance degradations and improvements
To compare the performance of the different variations I have run with perf-math repo from @seantalts $10^4$ evaluations of the lgamma function with 10 replications. The compared cases are

old_boost: boost lgamma with the old policy which does propagation from double to long double
boost: boost lgamma with the new policy avoid the over propagation
std: the old std::lgamma which is not reentrant safe, but currently in develop
cstd: the lgamma_r which is used on macOS and Linux which is reentrant safe

Note that these tests use the patched boost lgamma already.
As the boost implementation works differently for small (<1E25) and large (>1E25) arguments, I did two variations of the test.
# small arguments
BM_old_boost_median     843457 ns     842196 ns        732
BM_boost_median         266367 ns     266276 ns       2667
BM_std_median           204035 ns     203901 ns       3324
BM_cstd_median          189103 ns     189051 ns       3583

# large arguments
BM_old_boost_median     888854 ns     886906 ns        736
BM_boost_median         135407 ns     135268 ns       5415
BM_std_median           207349 ns     207202 ns       3379
BM_cstd_median          204118 ns     203722 ns       3449

Thus, on Windows we will have a ~30% performance degradation whenever no STAN_THREADS is set since we could use std::lgamma if we wanted to branch on Windows between thread use or not. However, for very large arguments the boost lgamma is even significantly faster than any other implementation.
On macOS and Linux we are getting a 5% performance increase, since the lgamma_r is a little faster than the current std::lgamma from these tests.
Checklist

[X] Math issue #1250

[X] Copyright holder: Sebastian Weber
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[X] the basic tests are passing

unit tests pass (to run, use: ./runTests.py test/unit)
header checks pass, (make test-headers)
docs build, (make doxygen)
code passes the built in C++ standards checks (make cpplint)

[X] the code is written in idiomatic C++ and changes are documented in the doxygen

[X] the new changes are tested

I don't see why using Boost's version unconditionally is warranted. If compilers are generally exhibiting poor performance when std::lgamma is called in parallel, then we can use Boost's version conditional on STAN_THREADS being defined. If it is just the Mac compiler, then we can also condition on that. Why would we want to take a performance hit and extra compilation time in the non-parallel case?

@bgoodri:  The issue says "I changed the default policy of boost to not do this propagation and now the speed is much closer to the std version."  So I don't think speed's an issue.
The motivation for not ifdef-ing all over the place is just code maintainabilty (which includes readability).  It wouldn't be difficult to add the ifdefs in this case, but there's always an ongoing maintenance cost in that the code gets more complex and harder to read.  That can be justified for big performance gains.

@serban-nicusor-toptal Could you please have a look at what is going on here? The tests are failing because files which are in boost are not found. These files are certainly in my local boost copy part of stan-math so that I am a bit confused.
Thanks.
@bgoodri As I have written... the speed is now much closer to the std version. From how I read to docs of boost for this, then the precision of boost is better than the std version - so if we want to match the speed we could limit the number of digits which can be assumed to be correct in order to get even more speed. Moreover, this boost policy is not systematically invoked for all of our boost math calls - gamma_p and gamma_q  do not use the policy, for example. Thus, we can even presumably speed things up for those functions by switching them to use the policy (those policies can be turned on using compiler flags as well if we want to find out without too much hassle).
And after all, we will see the performance benchmark report once the test pipeline finishes. Then we will see what effect this change has.

@serban-nicusor-toptal Mind looking into it this one:

The Jenkins link here in github points to the wrong place
More importantly, the Windows tests failed for reasons which look like system setup problems

Thanks!

@wds15 Hey, I've enabled windows ec2 agents and they may miss some dependency, will fix a bit later.
Thanks for pointing out the url thing, I will fix that asap.

Bump. I just had another Windows hickup. Is this stable right now, should I keep hammering things or just wait until this is settled?

Everything should be fine now.
Can you tell me what happened now ?

@serban-nicusor-toptal the second windows hiccup looked like the first one - maybe the logs have more details, but it seems as if the trouble as gone away.

Sounds good. If anything else occurs please keep me updated!
Thanks

The precision of the lgamma from boost should be a bit better than the std version, but not much as I have configured the boost version to use double precision only (instead of going to long double precision inernally which is the default).
The performance tests posted above do test the performance, no? While I don't have an overview which model uses specifically lgamma - that function is used in a few models.
So what is wrong with the above test?
(also... we should probably accept a bit of performance loss if it exists, I think - given that the alternative is to pour ifdefs into out codebase and make the use of boost vs std dependent on use of STAN_THREADS or not).

Let me take another look

@wds15

that function is used in a few models.

sure, then I'm not worried about sampling precision

The performance tests posted above do test the performance, no?

Yeah this is correct. Which one of these performance test are using the lgamma function? I'll go ahead and approve this, as you said:

given that the alternative is to pour ifdefs into out codebase

Also - can we contribute files to the stat_comp_benchmark and have them tested? Would be cool to test more than what's up there (glm regression, etc).

@drezap Thanks for reviewing this. You should bring up the request for more or different models on the forum, I think. We need to balance test coverage numerical cost vs testing needs as usual. Ideally an optional more extended set can be tested, but that is work to implement that.
@syclik Are u fine with merging? You were pretty strong about adding a test which fails for the std lgamma. To avoid in the future the performance regression wrt to non-scalability of the std version, I chose to add documentation. Adding a test for the scalability would require to test if on 4 cores the speed of a reasonably sized problem is about 4x as fast. The reason to avoid this is

this test may only hold on macOS where I have seen this problem - the exact behaviour wrt to threading depends on the OS and the compiler as it is an implementation detail according to the standard. Given that things depend on the implementation we cannot be sure that the test fails on all platforms whenever we want it to fail (a Linux developer may potentially not notice at all depending on the Linux implementation of lgamma).

running such a test requires that 4 cores must be idle when running the test. If that is not the case, then the test will likely fail since no 4x speedup will be measured in that case. Thus, the testing codebase gets unstable to some degree as it depends on the computing environment of the machine.

On Jun 19, 2019, at 6:16 AM, wds15<EMAIL_ADDRESS>wrote:
The precision of the lgamma from boost should be a bit better than the std version, but not much as I have configured the boost version to use double precision only (instead of going to long double precision inernally which is the default).

I think a lot of the Boost functions may be doing that implicitly.  I'm wondering if there are other things we can speed up in the same way (not in this issue, of course).

The performance tests posted above do test the performance, no? While I don't have an overview which model uses specifically lgamma - that function is used in a few models.
So what is wrong with the above test?
(also... we should probably accept a bit of performance loss if it exists, I think - given that the alternative is to pour ifdefs into out codebase and make the use of boost vs std dependent on use of STAN_THREADS or not).

How much performance loss is there?  It would be ideal to test the function out of the context of an end-to-end model where all sorts of things are going on.  It's not autodiff, either, just the base function.  So a really easy standalone test is just to add up a bunch of outputs applied to random inputs in a loop. Don't time the random number generation.
std::vector random_inputs = ...;
start timer
double y = sum(lgamma(random_inputs));
end timer
print output (necessary to prevent sum from being compiled away) and run time

Just a heads up that this PR needs a merge with fresh develop due to the meta refactor in #1268

@syclik & @bob-carpenter So it looks like we will pay a 22% performance degradation for lgamma here if we merge this in the current form (on macOS). What do you suggest:

live with it and merge
evaluate if lowering the precision targets of boost brings us somewhere
pour some ifdefs in here to use std for non-threading code and boost for threading code (those ifdefs will pop-up in some of the tests as well!)

Not sure what is best (I would probably do 1 and live with it given the other options are not nice either).

@bob-carpenter Since you asked if we have more possibilites to speed up boost math functions by avoiding propagation from double to long double internally. I did look into our codebase and it did not take long until I found boost math functions which do not respect the boost_policy_t. Those functions will thus continue to propagate to long double. We can either catch all of these instance or just slab on the compiler flag -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false during compilation. That will globally disable the propagation, see here.

Ok... going down to 9-10 digits accuracy (base 10) is what is needed to make boost lgamma match the performance of the std version...does not sound convincing to me.

@syclik Are u fine with my suggestions here? I updated the PR comments to reflect the changed solution. Once you are ok with using lgamma_r I will look into how to make test-math-dependencies to reject changes here.
@bgoodri Any thoughts from your end? Can you manage an additional define when compiling Stan programs in R?

We can always have more defines. What happens on icpc?
On Thu, Jun 20, 2019 at 4:37 PM wds15<EMAIL_ADDRESS>wrote:

@syclik https://github.com/syclik Are u fine with my suggestions here?
I updated the PR comments to reflect the changed solution. Once you are ok
with using lgamma_r I will look into how to make test-math-dependencies
to reject changes here.
@bgoodri https://github.com/bgoodri Any thoughts from your end? Can you
manage an additional define when compiling Stan programs in R?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAZ2XKR4ITCIPHLHBKKQ7L3P3PS7TA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODYGRJIY#issuecomment-504173731,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAZ2XKVNTICDYYLRLM5KZYDP3PS7TANCNFSM4HRGN7BA
.

icpc supports all that GNU specific stuff and I just tried it - works for me with intel compiler 15.0.3
Hopefully I can trace down what is happening on windows mingw where there have been problems in the last round.

Thanks much for evaluating.
You mean it takes 22% longer with the new code?  Is the
accuracy the same?
If so, I don't think that's a huge deal as lgamma is
rarely a bottleneck in models.  It might slow down
some hierarchical gamma or beta models.
I wouldn't trust ifdefs on this unless we can condition on
some environmental variable that guarantees single-thread
execution.

On Jun 20, 2019, at 3:23 PM, wds15<EMAIL_ADDRESS>wrote:
@syclik & @bob-carpenter So it looks like we will pay a 22% performance degradation for lgamma here if we merge this in the current form (on macOS). What do you suggest:
• live with it and merge
• evaluate if lowering the precision targets of boost brings us somewhere
• pour some ifdefs in here to use std for non-threading code and boost for threading code (those ifdefs will pop-up in some of the tests as well!)
Not sure what is best (I would probably do 1 and live with it given the other options are not nice either).
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

Thanks for diving into this.  This seems like a reasonable option if we can live with reduced precision.  The main places we need precision are long sequences of operations.  We don't need a huge amount of precision in any single term that's added to the log density target.
There are a lot of Boost functions we call without adding policy arguments, too.  I found some recently in the autodiff testing.

On Jun 20, 2019, at 3:31 PM, wds15<EMAIL_ADDRESS>wrote:
@bob-carpenter Since you asked if we have more possibilites to speed up boost math functions by avoiding propagation from double to long double internally. I did look into our codebase and it did not take long until I found boost math functions which do not respect the boost_policy_t. Those functions will thus continue to propagate to long double. We can either catch all of these instance or just slab on the compiler flag -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false during compilation. That will globally disable the propagation, see here.

How much accuracy is in the std:: version?  How much is even possible (within reason, of course) for something like lgamma?

On Jun 20, 2019, at 3:55 PM, wds15<EMAIL_ADDRESS>wrote:
Ok... going down to 9-10 digits accuracy (base 10) is what is needed to make boost lgamma match the performance of the std version...does not sound convincing to me.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

Those precision questions are really hard and very laborious to answer... I don't know. My impression is that boost is better controlled in that regard.
Sorry for the slight mess here, but I had been exploring this a bit.
To my understanding now we could use on POSIX systems the lgamma_r function which is even 4% faster than the current std::lgamma which we have - that's great. BUT: the mingw gcc compiler on windows does not have this POSIX support. So on windows I would have to stick to the boost version which is 22% slower than the std one.
This way we would always use a thread safe reentrant version in any case - and on POSIX systems (Linux & Mac) we get a speedup of 4%. Is that a plan?

On Jun 21, 2019, at 1:11 PM, wds15<EMAIL_ADDRESS>wrote:
Those precision questions are really hard and very laborious to answer... I don't know. My impression is that boost is better controlled in that regard.

I know.  I thought you figured out how to do it since you reported precision in another message.  Was that just because Boost let you set the precision or because you could compare a few values vs. the high-precision setting?

Sorry for the slight mess here, but I had been exploring this a bit.

It's good to get this stuff out in the open and see what's going on.
Not your fault lgamma's an odd edge case that just happens to be very
important for normalizing constants in stats.

To my understanding now we could use on POSIX systems the lgamma_r function which is even 4% faster than the current std::lgamma which we have - that's great. BUT: the mingw gcc compiler on windows does not have this POSIX support. So on windows I would have to stick to the boost version which is 22% slower than the std one.
This way we would always use a thread safe reentrant version in any case - and on POSIX systems (Linux & Mac) we get a speedup of 4%. Is that a plan?

I wouldn't want to add special code to get a 4% speedup in lgamma on Linux/Mac.

The boost doc has some info on the precision which is where I took this from.
So the POSIX C-implementation of lgamma_r is ~4% faster than std::lgamma_r and the boost::math::lgamma is 22% slower than the std version.
Comparing the POSIX C-implementation lgamma_r with boost shows that the boost version is 34% slower. This is quite significant and you will notice that end-to-end if you do Poisson regression with normalization. So I think we should take advantage of lgamma_r on POSIX systems.
For details see the timings posted above in the PR description.

The unsupported part of Eigen has a lgamma_r implementation that you could
use on Windows.
On Fri, Jun 21, 2019 at 2:53 PM wds15<EMAIL_ADDRESS>wrote:

The boost doc has some info on the precision which is where I took this
from.
So the POSIX C-implementation of lgamma_r is ~4% faster than std::lgamma_r
and the boost::math::lgamma is 22% slower than the std version.
Comparing the POSIX C-implementation lgamma_r with boost shows that the
boost version is 34% slower. This is quite significant and you will notice
that end-to-end if you do Poisson regression with normalization. So I think
we should take advantage of lgamma_r on POSIX systems.
For details see the timings posted above in the PR description.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAZ2XKQIGYMDV675MJSSTBDP3UPRLA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODYJJNZA#issuecomment-504534756,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAZ2XKXACBKSBLTJEX7U4XDP3UPRLANCNFSM4HRGN7BA
.

@bgoodri I haven't found that Eigen implementation. Not sure which one you mean.

This is ready for review @syclik ?

@serban-nicusor-toptal Could you please have a look at the windows configuration? I am getting on a Windows box
echo CXX=clang++-6.0 -Werror > make/local
which will just not work since we use g++.
To avoid it I modified the Jenkinsfile to not use the env.CXX variable - but that variable seems to contain the wrong thing.
You can see the problem in the last failed build 17 of this PR.
Thanks.

@serban-nicusor-toptal Ok, now I am out of ideas. The windows machine appears broken to me. It would be great if you could check.

Hey, looking into it now.
Spawning a machine to see if the var is wrong in the template.

Is there anything wrong that I should change from there ?

Yes... we use g++ and not clang++-6.0 as CXX.
This was not a problem so far, because the windows tests just used the default compiler on windows which is g++.
What is confusing is that only the Windows threading tests fail, but not the other one - I did not see a difference in the Jenkinsfile though.
Thanks for checking this out.

Okay, will update the ami template and rebuild this job.
On Thu, 04 Jul 2019, 7:21 PM wds15<EMAIL_ADDRESS>wrote:

Yes... we use g++ and not clang++-6.0 as CXX.
This was not a problem so far, because the windows tests just used the
default compiler on windows which is g++.
What is confusing is that only the Windows threading tests fail, but not
the other one - I did not see a difference in the Jenkinsfile though.
Thanks for checking this out.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=ALR77DCNBHLFARHCGEJKHN3P5YWTNA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODZH4CYQ#issuecomment-508543330,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ALR77DCTWMBLGGUXRQOMND3P5YWTNANCNFSM4HRGN7BA
.

@serban-nicusor-toptal updating the ami template did not solve it. Same errors.

The only "error" found in the logs is this:
Running on EC2 (Windows) - Build, Docker, WSL (i-0cea8c7f6ebab4b1a) in J:\Jenkins\workspace\Math_Pipeline_PR-1255

J:\Jenkins\workspace\Math_Pipeline_PR-1255>attrib -r -s /s /d
File not found - J:\Jenkins\workspace\Math_Pipeline_PR-1255\*.*

A test is ran like the following:
g++ -std=c++1y -m64 -Wall -Wno-unused-function -Wno-uninitialized -Wno-unused-but-set-variable -Wno-unused-variable -Wno-sign-compare -Wno-unused-local-typedefs     -O3  -I . -I lib/eigen_3.3.3 -I lib/boost_1.69.0 -I lib/sundials_4.1.0/include -I lib/gtest_1.8.1/include -I lib/gtest_1.8.1   -D_USE_MATH_DEFINES  -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION    -DGTEST_USE_OWN_TR1_TUPLE  -c -MT test/unit/math/fwd/scal/fun/atan2_test.o -M -E -MG -MP -MF test/unit/math/fwd/scal/fun/atan2_test.d test/unit/math/fwd/scal/fun/atan2_test.cpp

No test report files were found. Configuration error?

Besides this there aren't any errors that may point to a reason for this failure.

@syclik This is ready now for review.

bump... it would be awesome to get this in time for 2.20 - this PR speedups up models which use the log gamma in a noticeable way (should be quite a few) and we get a reentrant safe version along with scalable behavior when threading is turned on.
In case there is a need to iterate (which I hope is not), then I would try to have short turn-around.

Thanks for your comments.
The -D_REENTRANT thing is more or less a standard thing as I understood. That should be safe to use on Linux & macOS... but since there are apparently exceptions, I wanted to leave the user of Stan-math a choice. So if for whatever reason -D_REENTRANT does not bring up a lgamma_r, then the user can drop this definition and stan math will go back to the always available boost version.
The lgamma_r is something you really want to use due to speed advantages, but just in case things are odd we should leave an option to go with the boost thing (noticeably slower though... and with the long double precision way too slow).
Does that make sense?

@syclik I forgot: I also do not think we can assume that someone uses stan-math with our makefiles - or can we? So whenever someone is not using our makefiles then the _REENTRANT thing is not necessarily defined such that your suggestion to always define it and then assume that lgamma_r is there would fail. We would force everyone to enable -D_REENTRANT.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Given that this PR changes the default policy of boost to avoid promotion of double to long double internally this has side effects on all functions using it. These are
./prim/scal/fun/rising_factorial.hpp
./prim/scal/fun/falling_factorial.hpp
./prim/scal/fun/digamma.hpp

The digamma is probably the key one here as its more often used. The digamma is involved in the CDF of the negative binomial. To evaluate the impact of the change I did run an example model with and without the changed policy. The result is that the Stan inference is basically unchanged, but the program runs twice (!!) as fast. See below:
data {
  int X;
  int x[X];
}
transformed data{
  int U[X];
  for(i in 1:X) U[i] = 150;
}
parameters{
  real mu;
  real<lower=0> sigma;}
model{
  real lp[2];

  lp[1] = neg_binomial_2_lpmf(x | exp(mu), sigma);
  // calls internally the regularized incomplete beta function which
  // we compute with boost math (digamma)
  lp[2] = - neg_binomial_2_lcdf(U | exp(mu), sigma);

  target += sum(lp);
}

Data file (rename to .R file please)
negbin_data_R.txt
The run with the released 2.20 version gives me this summary:
vanilla (current default policy with double to long double propagation)
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (57) seconds, 57 seconds total
Sampling took (54) seconds, 54 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.1e+00  -4266  -4263  -4263    424      7.9  1.0e+00
accept_stat__    0.93  2.6e-03  8.9e-02   0.74   0.96    1.0   1185       22  1.0e+00
stepsize__       0.64  7.8e-16  7.8e-16   0.64   0.64   0.64    1.0    0.019  1.0e+00
treedepth__       2.1  1.7e-02  5.5e-01    1.0    2.0    3.0   1027       19  1.0e+00
n_leapfrog__      4.9  7.7e-02  2.4e+00    1.0    3.0    7.0   1011       19  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500      9.3      nan
energy__         4265  7.9e-02  1.6e+00   4263   4264   4268    402      7.5  1.0e+00
mu                4.6  6.9e-04  1.5e-02    4.6    4.6    4.6    487      9.0  1.0e+00
sigma              11  3.3e-02  7.5e-01    9.3     11     12    521      9.7  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

fast (same program, but adding CXXFLAGS+=-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false to make/local which does the same override and avoids the propagation):
Inference for Stan model: negbin_trunc_model
1 chains: each with iter=(1000); warmup=(0); thin=(1); 1000 iterations saved.

Warmup took (28) seconds, 28 seconds total
Sampling took (22) seconds, 22 seconds total

                 Mean     MCSE   StdDev     5%    50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4264  5.1e-02  1.2e+00  -4266  -4263  -4263    532       24  1.0e+00
accept_stat__    0.87  4.8e-03  1.7e-01   0.51   0.93    1.0   1169       52  1.0e+00
stepsize__       0.85  3.1e-15  3.1e-15   0.85   0.85   0.85    1.0    0.045  1.0e+00
treedepth__       2.0  2.0e-02  6.0e-01    1.0    2.0    3.0    894       40  1.0e+00
n_leapfrog__      3.7  7.1e-02  2.0e+00    1.0    3.0    7.0    795       35  1.0e+00
divergent__      0.00  0.0e+00  0.0e+00   0.00   0.00   0.00    500       22      nan
energy__         4265  7.3e-02  1.6e+00   4263   4264   4268    477       21  1.0e+00
mu                4.6  7.2e-04  1.6e-02    4.6    4.6    4.6    509       23  1.0e+00
sigma              10  3.2e-02  7.5e-01    9.2     10     12    546       24  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Note the very pleasant doubling in  speed, but results are really the same to my eyes.

Ups...there was a GitHub system problem and now we have my last comment 5 times posted.

@serban-nicusor-toptal There seem to be again windows hiccups. Now I see during windows testing (already for the 2nd time):
J:\Jenkins\workspace\Math Pipeline\PR-1255: Unable to determine if root directory exists

We should really think about cutting back on Windows testing - it is a major source of problems.

Steve and Serban gave me access to a Windows AWS instance and I ran all the commands Jenkins runs
attrib -r -s /s /d
echo CXX=g++ -Werror  1>make/local
echo CXXFLAGS+=-DSTAN_THREADS  1>>make/local
runTests.py -j16 test/unit -f thread --make-only
runTests.py -j16 test/unit -f thread
runTests.py -j16 test/unit -f map_rect --make-only
runTests.py -j16 test/unit -f map_rect

and had no issues running them. Everything passed on the exact same instance as the Jenkins tests are being run on.
The last idea I had was that the spaces in the "J:\Jenkins\workspace\Math Pipeline\PR-1255:" are causing issue, but the Stan Windows tests fail with the same errors and dont have spaces in them so I guess not :/.
The Gelman Windows machine at Columbia that is used if AWS Windows instances are turned off, is having issues with the Permission error, which is I guess a sign of no disk space.
I agree that this is really frustrating... Can we have someone look at Gelman Windows machine so we at least have a backup until we find the cause fo this "unable to find root folder" error?

I think @serban-nicusor-toptal must be experimenting with the Windows EC2 fleet again. Nic, would you mind turning off the windows label from those machines and creating a new Math test job that uses an alternative label so that there's no interference with pull requests and the "production" continuous integration suite? Thank you.

Turning off AWS instance probably wont help since the Gelman Windows machine is not working either.

Hey, J:\Jenkins\workspace\Math Pipeline\PR-1255: Unable to determine if root directory exists is an aws ec2 instance issue, while we we're testing something messed up probably, I will fix it as soon as I get home.
Gelman machine does not have any issues, I've checked it last night.
@seantalts Disabled to not interrupt the master/prs. Will create a separate job.
@rok-cesnovar @SteveBronder thanks again for the help with this !

Thanks @serban-nicusor-toptal and @rok-cesnovar for your work on this! Windows is such a pain and bottlenecks a bunch of our builds, so hopefully we can figure something out.

I think the gelman machine is reportin Permissiom denied errors.

I don't remember how many times I already restarted, but I am getting once again on Windows
J:\Jenkins\workspace\Math Pipeline\PR-1255: Unable to determine if root directory exists

It seems to be working now, 4 green light on Always run tests part 2.

@syclik This is good to review now.
Thanks for everyone handholding the Jenkins. Yack.

Huh?? Yesterday this PR was all green and now the Jenkins runs again without any activity?

(stat_comp_benchmarks/benchmarks/arK/arK.stan, 1.26)
(stat_comp_benchmarks/benchmarks/garch/garch.stan, 1.11)

Are the tests always run against the current develop version?
I see these in my other PRs as well, ones where we don't really touch the things behind either of these models.
@serban-nicusor-toptal can you link me to the file that run the tests? I had a hard time finding it

I think that those tests are just too fast and we are seeing a lot of noise here. At least this is what I recall. So maybe we drop these from being tested?

If they are v simple maybe we should just run them more than once

@syclik are you reviewing this or can I have a look? Dont want to duplicate if you are on it.

@rok-cesnovar, thank you for asking. Please feel free to review anything if there's a delay from me. I'm sure we'll all reach out for help where it's uncertain.
I know the binomial coefficient log function gets used a lot and we've had numeric stability issues there before. I'm worried that these side effects will prevent users from fitting models that they are able to fit now because once you're in those tails, you can never escape due to overflow.

If that’s the case, why didn’t you add a new test?
That would have made it much clearer to a reviewer and you wouldn’t have
needed to add a lengthy comment in the code. This sort of comment adds
additional work in the future: why was it so important it needed to be
called out? That’s why I had interpreted it as being an undocumented side
effect now.
I think the quick fix is to add back the test you removed (unless it has no
purpose, although it sounds like it should be testing an overflow
condition) and remove your comment. What do you think?
On Thu, Aug 1, 2019 at 5:24 AM wds15<EMAIL_ADDRESS>wrote:

@wds15 commented on this pull request.
In test/unit/math/prim/scal/fun/binomial_coefficient_log_test.cpp
https://github.com/stan-dev/math/pull/1255#discussion_r309603878:

@@ -27,7 +28,14 @@ TEST(MathFunctions, binomial_coefficient_log) {
}

test_binom_coefficient(1e9, 1e5);

test_binom_coefficient(1e50, 1e45);

// NOTE(SW): 2019-05-31 replaced large number test as

The large number test did always overflow (so it still does on develop).
The std lgamma returned Inf while the boost returns NaN. I need to check
the lgamma_r, but I would guess it returns Inf (the lgamma_r only differs
in how the sign is stored of the argument). I changed the test so that we
avoid the overflow and actually do a meaningful test here.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AADH6F3XLDBCTSDLCRG7PI3QCKTTHA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFWFIHK3DMKJSXC5LFON2FEZLWNFSXPKTDN5WW2ZLOORPWSZGOCAIDNUI#discussion_r309603878,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AADH6FYCZIHRSQHE5FASN4TQCKTTHANCNFSM4HRGN7BA
.

We cannot add back the test, because Inf comparisons seem to work for google test, but nan comparisons don’t work.
I was very confused by the test myself. From the history the test was introduces by bob a while ago. I guess he wanted to test the large number branch of the binomial coefficients function. I guess it wasn’t noticed that the test did overflow and never tested actual numbers.
I am happy to change the comment to what we think makes things easiest. I have put the comment there to not forget...so it’s good that we discuss it (but I just forgot which is why I put comments here and there).

@wds15, I think it’s clear that you’re introducing behavior that was not intended and it will have side effects that we don’t even know the extent of. Is there any way we can have the benefits you are looking for without changing existing behavior everywhere else?
NaN and Inf are very different, especially as it relates to the autodiff stack. If there’s a test that shows the function is expected to return infinity, then we should do that. If that behavior needs to change, please make that a separate issue that we can discuss independently of this issue.

I donˋt think it makes sense to insist on tests as they ar just for the sake of it. If we change them, then there must be very good reasons, of course. The test in question about the binomial coefficient was testing Inf behavior on accident..maybe @bob-carpenter could maybe comment here if he wanted to test inf behavior or just large numbers (now we test large numbers), but I really doubt that inf testing was intended.
The culprit of this is not defined behavior of lgamma for too large values. In the past we moved from boost lgamma to std. This change at the time reduced the precision of lgamma and it changed the behavior when overflowing (maybe even unnoticed, so this pr is far more transparent, just btw). I am not aware of a convention in stan math for this. Moreover, if we insist on inf being returned, then we have a problem with the other boost functions like the digamma.
So, if we want to test for overflow behavior of lgamma, then we should test that explicitly and not via the binomial coefficient test...but then the other boost functions must be checked as well.
If you want the same overflow behavior, then we have to find a large number like 1e20 where the lgamma implementations are about to flow, but we catch that with a bloody if statement and then return inf.
Finally, could you please explain why nan or inf makes such a difference? Both will lead to a rejection of the draw in stan, no? The arithmetic with inf is super fragile from what i have seen.

Is this for the performance regression or longer term?
When I recently updated the autodiff test framework, I added extensive tests for all of the vectorized unary functions that include boundary tests for +inf, -inf, and NaN.  It makes sure primitive (double) version produces the same value and and finite-difference derivatives as autodiff results up to third order (specifically: var, fvar<double>, fvar<fvar<double>>, fvar<var>, fvar<fvar<var>>).  Some boundaries are not checked because we often provide derivatives if the limits exist from one side even though finite differences will blow up.  We may want to refine all of that boundary condition behavior going forward.
I think @syclik needs to decide on policy here.
Also, could someone please review my pull request.  It's been open for 2 months and is going to get stale.  The intent was to define behavior and check all boundary conditions like these automatically.

@wds15, I took time to respond to your reply because I wanted to make sure I responded without haste.
The issue and PR states that this fixed lgamma and does not mention binomial_coefficient. From a review point of view: this PR should not be accepted because it has unintended side effects. As a general principle, fixes to one part of the code shouldn’t alter other parts of the code. That’s why we have tests.
This is not suggesting that the code is frozen. If you read my response, I suggested adding a new issue if you’ve discovered a new problem. (In this case, I think we would reject what you’re proposing for change, but we should evaluate that separately.)
@bob-carpenter suggested a policy. I believe he meant for boundary conditions. We can and should discuss that separately.
For PRs, we should not change behavior for other parts of the code without deliberately doing it and without deliberation. Just imagine if this was someone else’s PR. You really wouldn’t want someone changing a test to match behavior to fix something unrelated and just let that slide. That’s not a good PR policy to have.
So... if you need to have a conditional to guarantee the behavior stays elsewhere, then use a conditional. If you need to change the behavior, then file a new issue that we (the community of Math developers) can discuss. If you need to close this PR to rewrite the implementation and come up with a clearer top-level description, that’s ok too. If this were a bug, I’d still expect a new issue to outline the bug.
For specifying the PR and the issue, please try to write side effects down. If I had read this in the description, I would have told you that this isn’t something that should be changed lightly and could have saved us a lot of time. (I had to dig to find a change that was snuck in the PR... that’s not good practice.)
Those are a lot of options. What would help you get this PR up to an implementation that could be merged into the codebase? If I were to suggest a path, I’d say start with fixing lgamma while maintaining the behavior that existed everywhere else. That’s the quickest path to merging something.

Let me amend that general principle to say:

Fixes to one part of the code shouldn't break other parts of the code.

Almost anything we change at the math library level is going to have repercussions up the function chain.  For instance, any change to lgamma is going to percolate through to binomial_coefficient, beta_lpdf, etc.
That's usually a good thing when we fix bugs.
The issue here is presumably that we are changing behavior of lgamma, which has a lot of repercussions. lgamma is used on 77 lines based on
> grep -R "lgamma(" stan/math | wc

Given that there's a bug due to the interaction of multithreading and the existing lgamma implementation, how do you want to proceed?  I believe this PR is the attempt to fix lgamma.  Do you have other ideas for how to fix it?
P.S.  I'd like to stay away from language like "snuck in" as it implies an intention to deceive that I doubt was there in this PR.

The culprit is really that we have not a well defined behavior of lgamma for overflowing. We only test overflowing behavior implicitly via the binomial coefficient. The binomial coefficient test works if overflowing of lgamma returns inf...which it doesn’t do if we use boost on windows. The boost variant returns nan when it overflows, whih breaks the binomial coefficient test. As we used boost earlier I thought things are fine and I just change the binomial coefficient test to actually test large numbers and not inf.
So i suggest to do:

enable the binomial coefficient test for inf of lgamma again
explicitly test lgamma to overflow with inf for huge arguments ...which is the correct way for me to ensure this and not implicitly (but we can leave the implicit test due to policy)
either use a conditional in lgamma for the overflowing to be right or configure boost to return inf on overflow (i looked into that already and failed last time, will look again)

(I was thinking changing lgamma is easy, but i am again trapped in a 90+ Comments pr, sigh)
I still don’t understand why inf vs nan is such a big deal and I would very much like to understand...I just don’t know (nan or inf log-lik or gradients lead to rejection to my understanding)

On Aug 8, 2019, at 4:43 PM, wds15<EMAIL_ADDRESS>wrote:
The culprit is really that we have not a well defined behavior of lgamma for overflowing. We only test overflowing behavior implicitly via the binomial coefficient. The binomial coefficient test works if overflowing of lgamma returns inf...which it doesn’t do if we use boost on windows. The boost variant returns nan when it overflows, whih breaks the binomial coefficient test. As we used boost earlier I thought things are fine and I just change the binomial coefficient test to actually test large numbers and not inf.
So i suggest to do:
• enable the binomial coefficient test for inf of lgamma again
• explicitly test lgamma to overflow with inf for huge arguments ...which is the correct way for me to ensure this and not implicitly (but we can leave the implicit test due to policy)
• either use a conditional in lgamma for the overflowing to be right or configure boost to return inf on overflow (i looked into that already and failed last time, will look again)
(I was thinking changing lgamma is easy, but i am again trapped in a 90+ Comments pr, sigh)
I still don’t understand why inf vs nan is such a big deal and I would very much like to understand...I just don’t know (nan or inf log-lik or gradients lead to rejection to my understanding)

We just want to get the limits right, o that lgamma(inf) should be inf, for example.  I don't think it's a huge deal.
A condition to return infinity is fine.
In practice in Stan, it's not such a big deal as it's going to wind up causing numerical problems and rejections either way in most cases.  That is, our autodiff arithmetic won't survive a round trip exp(log(0)) even though log(0) = -inf and exp(-inf) = 0.

@wds15: yes, I think what you outlined is good. Just find some way to implement lgamma without changing behavior to other parts of the code that are tested.
As I mentioned before, it's reasonable to change behavior, but that must come in a separate issue. To clarify, I don't believe you're intentionally trying to change behavior to do something malicious. (I never thought that way.) I do think you are changing behavior and not considering the consequences of this.
I will come up with an example of where NaN vs Inf matters, but it'll have to wait until after StanCon. We can discuss in person if you'd like. What @bob-carpenter mentioned is correct, but there are some tricky edge cases that I don't think you're considering. (If I recall, NaN in on the autodiff stack within the chain method will kill the whole computation even if that's not part of the gradient computation whereas Inf will not.)

Let's look at a slightly longer time horizon before considering what behavior is locked in.
The race condition was introduced in a commit on Dec 9, 2018 (don't know when it was merged), from @andrjohns.
https://github.com/stan-dev/math/commit/01bb6645815cd5546701c0575fcb4f38d4a41cca#diff-6d2e0c5e89ad16a152ce65335ca2e93c
with this diff:
- inline double lgamma(int x) { return boost::math::lgamma(x, boost_policy_t()); }
+ inline double lgamma(int x) { return std::lgamma(x); }

The current PR is merely returning to the behavior we had (other than MINGW32__ apparently) before the above commit with:
#if !__MINGW32__
  int sign = 1;
  return ::lgamma_r(x, &sign);
#else
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
#endif

P.S.  NaN can be made to matter if you use is_inf() or is_nan().  That'll come up in some boundary condition tests, where some functions might throw on NaN and let infinite values slide. Otherwise, it's going to cause the same kind of problems as sqrt(square(0)) [which introduces a 0 in the denominator in the derivative of sqrt(0)]. The method chain() is pure virtual; nothing dictates NaN behavior.  If you get a function called during chain() that checks boundaries, it might throw exceptions on either NaN or Inf.

Having lgamma(0) return Inf is correct mathematically in the usual
statistical context where its argument is assumed to be non-negative, what
the C++11 standard does, and is important, for example, with infinite sums
that have (l)gamma in their denominator and we want to terminate when
(l)gamma overflows to infinity.
On Mon, Aug 12, 2019 at 12:02 PM Bob Carpenter<EMAIL_ADDRESS>wrote:

Let's look at a slightly longer time horizon before considering what
behavior is locked in.
The race condition was introduced in a commit on Dec 9, 2018 (don't know
when it was merged), from @andrjohns https://github.com/andrjohns.
01bb664#diff-6d2e0c5e89ad16a152ce65335ca2e93c
https://github.com/stan-dev/math/commit/01bb6645815cd5546701c0575fcb4f38d4a41cca#diff-6d2e0c5e89ad16a152ce65335ca2e93c
with this diff:

inline double lgamma(int x) { return boost::math::lgamma(x, boost_policy_t()); }

inline double lgamma(int x) { return std::lgamma(x); }

The current PR is merely returning to the behavior we had (other than
MINGW32__ apparently) before the above commit with:
#if !MINGW32
int sign = 1;
return ::lgamma_r(x, &sign);
#else
if (unlikely(x == 0.0))
return std::numeric_limits::infinity();
return boost::math::lgamma(x, boost_policy_t());
#endif
P.S. NaN can be made to matter if you use is_inf() or is_nan(). That'll
come up in some boundary condition tests, where some functions might throw
on NaN and let infinite values slide. Otherwise, it's going to cause the
same kind of problems as sqrt(square(0)) [which introduces a 0 in the
denominator in the derivative of sqrt(0)]. The method chain() is pure
virtual; nothing dictates NaN behavior. If you get a function called
during chain() that checks boundaries, it might throw exceptions on
either NaN or Inf.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAZ2XKQKE6MUIBJIYBAMA73QEGCS3A5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4DAGPY#issuecomment-520487743,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAZ2XKXAMPHIG54YKQRYQUDQEGCS3ANCNFSM4HRGN7BA
.

Would it be sufficient to

add a condition to the double definition:

if (x == 0) return std::numeric_limits::infinity();

add a condition to chain() in the lgamma vari that if  value() == 0 then return what?  We need a derivative here.

add a condition to fvar implementation of lgamma that does the same thing as (2) for the derivative if value() == 0.

On Aug 12, 2019, at 12:26 PM, bgoodri<EMAIL_ADDRESS>wrote:
Having lgamma(0) return Inf is correct mathematically in the usual
statistical context where its argument is assumed to be non-negative, what
the C++11 standard does, and is important, for example, with infinite sums
that have (l)gamma in their denominator and we want to terminate when
(l)gamma overflows to infinity.
On Mon, Aug 12, 2019 at 12:02 PM Bob Carpenter<EMAIL_ADDRESS>wrote:

Let's look at a slightly longer time horizon before considering what
behavior is locked in.
The race condition was introduced in a commit on Dec 9, 2018 (don't know
when it was merged), from @andrjohns https://github.com/andrjohns.
01bb664#diff-6d2e0c5e89ad16a152ce65335ca2e93c
https://github.com/stan-dev/math/commit/01bb6645815cd5546701c0575fcb4f38d4a41cca#diff-6d2e0c5e89ad16a152ce65335ca2e93c
with this diff:

inline double lgamma(int x) { return boost::math::lgamma(x, boost_policy_t()); }

inline double lgamma(int x) { return std::lgamma(x); }

The current PR is merely returning to the behavior we had (other than
MINGW32__ apparently) before the above commit with:
#if !MINGW32
int sign = 1;
return ::lgamma_r(x, &sign);
#else
if (unlikely(x == 0.0))
return std::numeric_limits::infinity();
return boost::math::lgamma(x, boost_policy_t());
#endif
P.S. NaN can be made to matter if you use is_inf() or is_nan(). That'll
come up in some boundary condition tests, where some functions might throw
on NaN and let infinite values slide. Otherwise, it's going to cause the
same kind of problems as sqrt(square(0)) [which introduces a 0 in the
denominator in the derivative of sqrt(0)]. The method chain() is pure
virtual; nothing dictates NaN behavior. If you get a function called
during chain() that checks boundaries, it might throw exceptions on
either NaN or Inf.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAZ2XKQKE6MUIBJIYBAMA73QEGCS3A5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4DAGPY#issuecomment-520487743,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAZ2XKXAMPHIG54YKQRYQUDQEGCS3ANCNFSM4HRGN7BA
.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

@wds15: just FYI, this is what causes it to fail if we start returning NaN instead of inf: https://discourse.mc-stan.org/t/known-gradient-breaking-behaviours/9813/5 (I won't come up with a separate example unless you still need to see one.)
@bob-carpenter, the change was for issues that were clearly stated and had PRs approved for it: #1010 and #1055. At the time of the PRs was when the boundary conditions were tested.  We are trying to go in that direction with the Math library: be more consistent towards the std functions. There were no tests that would have failed for changing to std::lgamma at the time, so the implementation was changed and done so intentionally.
Re: the solution, I believe the first point was mentioned by @andrjohns in the original issue: https://github.com/stan-dev/math/issues/1010#issuecomment-445527864. For the derivative, are we going to have the same problem with digamma?

Is it just a matter of updating the PR text?
On Mon, Aug 12, 2019 at 21:07 Daniel Lee<EMAIL_ADDRESS>wrote:

@wds15 https://github.com/wds15: just FYI, this is what causes it to
fail if we start returning NaN instead of inf:
https://discourse.mc-stan.org/t/known-gradient-breaking-behaviours/9813/5
(I won't come up with a separate example unless you still need to see one.)
@bob-carpenter https://github.com/bob-carpenter, the change was for
issues that were clearly stated and had PRs approved for it: #1010
https://github.com/stan-dev/math/issues/1010 and #1055
https://github.com/stan-dev/math/issues/1055. At the time of the PRs
was when the boundary conditions were tested. We are trying to go in that
direction with the Math library: be more consistent towards the std
functions. There were no tests that would have failed for changing to
std::lgamma at the time, so the implementation was changed and done so
intentionally.
Re: the solution, I believe the first point was mentioned by @andrjohns
https://github.com/andrjohns in the original issue: #1010 (comment)
https://github.com/stan-dev/math/issues/1010#issuecomment-445527864.
For the derivative, are we going to have the same problem with digamma?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAGET3CESEC5FXWA7BTCCUDQEG7GXA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4DVQGQ#issuecomment-520575002,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAGET3BST6FWOABXMTZEKTLQEG7GXANCNFSM4HRGN7BA
.

Woha!!!
What a mess here: So the boost math lgamma turns out to be a really bad implementation, I think. I was proposing to avoid the internal double to long double propagation which happens in the default configuration of boost. Doing so speeds up the boost function up to an acceptable speed which is (still slower) comparable to the std::lgamma... BUT avoiding the propagation introduced overflowing of the boost lgamma at around 1E25 or so. This is due to a bad implementation of boost.... when you do allow for the propagation to long double, then the boost lgamma will be really slow, but it will return lgamma for up to 1E300 - just like the std::lgamma or ::lgamma_r whch also compute log gamma for up to 1E300 without problems!!! Getting results up to 1E300 makes defining overflowing behavior unneeded as 1E308 is the largest double there is.
Argh.
To get the best speed under all situations we could do this:

On Linux & macOS we use lgamma_r always (its fast and reentrant safe)... great. Done.
On Windows we use std::lgamma whenever STAN_THREADS is NOT defined while we use the boost::math::lgamma implementation when STAN_THREADS IS defined.

This would give the user always the fastest possible lgamma, but we have to juggle 3 calls here. In case we go with these 3 calls, then I would have to rethink the test-math-dependency check which does not work any more given this complexity.
My apologies for my wrong statements about the test_binomial_coefficient test for large arguments - that test does indeed test the correct thing.
(I will look into digamma of boost to see the state of affairs there)
If this is all confusing...it is! Maybe I quickly explain next week in person.

In case we go with these 3 calls, then I would have to rethink the test-math-dependency check which does not work any more given this complexity.

For this you have 2 options:

add an exception to the call_math_warnings function call to ignore the stan/math/prim/scal/fun/lgamma.hpp file in the check. This can be done using the fifth argument. I think this should do it:

@$(call math_warnings, 'thread-reentrant-safe', 'stan/math',\
             'std::lgamma',\
             'Replace std::lgamma with reentrant safe version boost::math::lgamma or lgamma_r on supporting platforms.',
             'stan/math/prim/scal/fun/lgamma.hpp')

dont add anything to test-math-dependencies in this PR but just make a post to #1294 explaining in short what we should check and what are the exceptions and I will add that once I do the rewrite. That should be ready the week following Stancon.

So I filed two issues for this for boost math. See here https://github.com/boostorg/math/issues/242 and here https://github.com/boostorg/math/issues/241 .  Let's see what boost has to say about this - they should really fix this, I think.
@syclik : The way I have coded the lgamma right now is this:
inline double lgamma(int x) {
#if !__MINGW32__
  int sign = 1;
  return ::lgamma_r(x, &sign);
#else
#ifdef STAN_THREADS
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
#else
  return std::lgamma(x);
#endif
#endif
}

NOTE: I was too tired yesterday to actually push the updates (forgot to commit), so the above is what I will push later today... sorry.
This way we branch for windows / no windows and we branch on windows for STAN_THREADS and not STAN_THREADS.
I know branching is not nice, but I don't see a better way. If we force to use boost lgamma on Windows regardless of STAN_THREADS, then we would make Windows users suffer from the slowness (given that the double to long double propagation has such nasty side effects I will take this out of this PR... fast is great, but overflowing is worse).
Let me know if that's good and then I update the PR description, etc..
@rok-cesnovar Thanks, I will then exclude the lgamma.hpp; that's probably the best thing to do.

good news: the boost folks appear to be very reactive... so maybe they fix matters quickly / hopefully. If boost gives us a fast and not overflowing version, then we should use that; let's see.

Looks like the boost folks have fixed the issue with the overflowing of lgamma whenever the double to long double propagation is turned off. See here:
https://github.com/boostorg/math/commit/fd07b121a619e0429e189edb4ff50f5b4e9efe58
Is it an option that I apply these few patches to our existing boost? If we do that then we can get a relatively fast lgamma from boost which we can use under Windows in all cases.
We could even consider to not take advantage of the additional speed from lgamma_r and instead use boost here as well (at the cost of ~10-20% performance loss; I would check that). The upside of doing so is that we

can use boost for all configurations and all settings => much easier to read code
the boost version gives us a platform-independent implementation of lgamma (whereas the std::* implementations are compiler + OS specific)

If we go with boost for most/all cases then we need to change the default policy of boost to avoid the propagation from double to long double. I have pointed out the performance loss in an issue to boost math developers (https://github.com/boostorg/math/issues/241) and they are even considering to change the default of boost math in an upcoming release. So doing this seems something quite sensible and (as I understand) the need to propagate from double to long double was justified long ago when matters where different (I can only guess as to why, but it sounds like this is a relict from the 32bit past, maybe).

Is it an option that I apply these few patches to our existing boost?

I would say that you can, since the precedent is that we already do that. Right now only this is patched : https://github.com/stan-dev/math/pull/1223/commits/b5f502a510e5239553aa8644c6a64cd378b8fc39 (explanation here) but it does set a "standard". There were more patches to Boost before but have since been reduced.
We just need to prevent any upgrades to Boost for version 1.70 and all other versions that would come before this patch is part of an official release.
We might actually be in luck there as Boost 1.71 is currently in beta and the official 1.71 release should come in a few days (time from beta to release for them is usually a month, its been a month and a half now since beta). Its just a matter if this patch makes it in. @wds15 would you mind asking them if this will get in?. If yes, we could also just wait for a week or two and upgrade to 1.71. If it wont make it in, we should prevent upgrades to 1.70 and 1.71 in the upgrade script https://github.com/stan-dev/math/blob/develop/lib/upgrade-boost.sh and that is it.

Traditionally we stayed with the boost version which is in the BH package on can. That is right now 1.69.0. If we patch boost we need to be careful, sure. So we should include a unit test which will trigger and overflow in case boost is not patched as a way to remind ourself of the need for the patch.

Ah yes, forgot about the BH & CRAN limitation, sorry. Adding the test you mentioned seems the cleanest and most effective solution.

We've patched Boost in the past. If it is just one file, it is no big deal.
On Wed, Aug 14, 2019, 4:12 AM Rok Češnovar<EMAIL_ADDRESS>wrote:

Ah yes, forgot about the BH & CRAN limitation, sorry. Adding the test you
mentioned seems the cleanest and most effective solution.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAZ2XKQ5HDD64ON2UQAS4KTQEO47BA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4IBMOI#issuecomment-521147961,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAZ2XKRUHLIU5RST7PKQSRTQEO47BANCNFSM4HRGN7BA
.

bump. This is good to review once tests turn green. PR description is updated to reflect all the nasty things one can bump into here.

@serban-nicusor-toptal could you have a look at the jenkins links, these don’t work right now. Thanks.

Uhm... ok... now it works...i was directed to a dashboard two h ago...not sure what was going on

Bumping this, are there any blockers at this point?

No. Everything is up to date as we discussed...just needs final review. This pr now patches boost to avoid overflowing at the lower precision.

Cool request who should review and so we can rock and roll!!

@syclik has an open "change requested" review... so it's on him? @rok-cesnovar did follow this closely as well, not sure if policy permits him to approve.
the main thing to do now, I think, is to verify that I back ported the boost lgamma patch correctly (I spend attention on to it).

I think Daniel should re-review this one, I think that makes sense. If he suspects he will not have time to do it within a reasonable time frame than he can assign me. I did follow the PR closely yes.
One thing I am wondering is what was the final number on the slowdown of boost::math::lgamma compared to the non-thread safe std::lgamma? So we can write down somewhere why we opted for boost:math::lgamma-only on windows and not do a  branch for the non-threaded case.
If you dont have easy access to a Windows machine I can run that, if anyone points me to a model & data to run.

The performance tests for the Windows question are basically done... I mean the numbers reported in the PR description are run on macOS, but they should round about hold for Windows as well.
So in Windows you will get BM_boost_median for the single-core performance which is somewhat (~30%) slower than the current BM_std_median for smallish numbers (1E25), but for large numbers the boost version is even faster. Maybe things are better when actually run on Windows, but I doubt that.
I already feel guilty of going to branching here, so I did not want to do another branching. Sure, one could use std::lgamma on Windows whenever STAN_THREADS is not set, but then I would also have to rethink the test-math-... thing and what if not.

@syclik If you could make a start separating out the side-issues (even just as a sketch) then that would be very helpful. I did try already to structure it nicely, but it failed apparently; so its most efficient if you can lay out the structure.
I can only run performance tests on Windows in a virtual machine; will take a moment. In case someone else can do that, the test is in https://github.com/wds15/perf-math under the branch lgamma and there the file lgamma.cpp.
We used boost math lgamma in the past. Now it is for sure a lot faster - as such I think its just fine. We are putting way to much burden on us developers to test every single configuration for optimal speed. This is an open-source project. No one gets paid for it. I find this too much effort. The mega burden on performance testing will not encourage future PRs to be reported with all their glory details (usually specific to the developers machine) as its easier to leave things in the dark. That's my opinion - more automated stuff is needed to warrant this level of coverage or we just accept that not every commit we make is an improvement in every aspect.

The results from my notebook running Windows 10 on an i7-5600U.
offset = 0.5 (is that considered small values performance, I honestly dont know)
Running lgamma.exe
Run on (4 X 2594 MHz CPU s)
CPU Caches:
  L1 Data 32K (x2)
  L1 Instruction 32K (x2)
  L2 Unified 262K (x2)
  L3 Unified 4194K (x1)
----------------------------------------------------
Benchmark             Time           CPU Iterations
----------------------------------------------------
BM_old_boost    1281636 ns    1260553 ns        533
BM_boost        1269738 ns    1255020 ns        498
BM_std           991667 ns     976562 ns        640

offset = 1E30
Running lgamma.exe
Run on (4 X 2594 MHz CPU s)
CPU Caches:
  L1 Data 32K (x2)
  L1 Instruction 32K (x2)
  L2 Unified 262K (x2)
  L3 Unified 4194K (x1)
----------------------------------------------------
Benchmark             Time           CPU Iterations
----------------------------------------------------
BM_old_boost    1261048 ns    1255020 ns        498
BM_boost        1236808 ns    1255580 ns        560
BM_std           968757 ns     976562 ns        560

Thanks. small and large is just below or above 1E25 in this context. So as predicted we see ~25% performance hit due to this change on Windows.

So to sum up we have the following options:

All options implement a thread-safe lgamma that we currently dont have but should have.

A (what is implemented in the PR now):
inline double lgamma(double x) {
#if !__MINGW32__
  int sign = 1;
  return ::lgamma_r(x, &sign);
#else
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
#endif
}

We get a 5% performance gain on Linux/MacOS and a 25% performance hit on Windows (see https://github.com/stan-dev/math/pull/1255#issuecomment-525627412).
B (the least "branchy" implementation - less maintanance):
inline double lgamma(double x) {
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
}

We get a 25% performance hit across the board. However we do have a 50% performance gain for very large arguments. See the top-level comment.
C (most "branchy" solution - hardest to maintain):
#if STAN_THREADS
#if !__MINGW32__
  int sign = 1;
  return ::lgamma_r(x, &sign);
#else
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
#endif
#else
return std::lgamma(x);
#endif

No performance hit, a 5% speedup on linux/macOS actually. But obviously difficult to read/maintain.
D (without lgamma_r):
#if STAN_THREADS
  if (unlikely(x == 0.0))
    return std::numeric_limits<double>::infinity();
  return boost::math::lgamma(x, boost_policy_t());
#else
return std::lgamma(x);
#endif

No performance hit, boost is a bit slower compared to lgamma_r on Linux/Mac (see top-level comment).
@syclik wrote on Discourse that we should follow this order of preference

Standard template library (STL)
Boost Math
CMath or other builtin C / C++ library that’s portable
custom implementation

STL is a no-go for the threaded case which means Boost is next in line in the order(option B above). If we put a huge weight on single-core performance we then have option D that is somewhere in between 1. and 2.
If we feel that lgamma_r performance shown in the top level comment is enough for CMath to overtake Boost in the order of preference then its either A (the one currently on this PR) or C, depending on where we stand in the "25% performance hit for Windows vs more maintenance" trade-off.
My personal view (does not carry much weight here, I am not the author or reviewer) is that if we follow the order of preference for implementations we should go with either D or B. I dont think the difference between boost and lgamma_r is enough to justify CMath prevailing over Boost.

C doesn't seem like it's that bad. sure there is some preprocessor stuff but it's not horrific and we can move back to a simpler version if the above becomes too much

For the record, I dont feel its "that bad" either. But if we have an order of preference for implementations we should stick to it or have a good reason to have a different order for a given function (for example if the threaded performance of lgamma_r is that much better than boost::math::lgamma then C might be the best).
Bottom line is that I think we have all the options for the final touch on this PR with the + and - for all of them.

There is one last consideration: The PR should also prevent future use of std::lgamma whenever threading is turned on... and this gives preference to any options which don't let any std::lgamma lurk into the codebase.

True, but we can exclude stan/math/prim/scal/fun/lgamma.hpp from the check to only allow std::lgamma be used in this file. If we really want to. But there is still some danger there off course.

@syclik If you could make a start separating out the side-issues (even just as a sketch) then that would be very helpful. I did try already to structure it nicely, but it failed apparently; so its most efficient if you can lay out the structure.

@wds15, I just added headings for each side effect. I think it could be clearer, but at least it breaks apart the different effects. One question: do any of these need to spawn a new issue?

We are putting way to much burden on us developers to test every single configuration for optimal speed. This is an open-source project. No one gets paid for it. I find this too much effort. The mega burden on performance testing will not encourage future PRs to be reported with all their glory details (usually specific to the developers machine) as its easier to leave things in the dark. That's my opinion - more automated stuff is needed to warrant this level of coverage or we just accept that not every commit we make is an improvement in every aspect.

I completely agree with you that we need to automate this. I thought I asked you for help in designing this, but if they wasn’t clear, please help. We should do this sooner rather than later.
Yes, everyone is a volunteer. But the burden isn’t solely yours; I’ve spent a lot of hours on this PR when the easy solution is to just say it’s not ready. And that’s just my time. Clearly, other developers are pitching in to help, much more than I am on this particular PR. I disagree that the burden is too high; if a developer wants to change the code base, the developer (and the community of developers) has a duty to make sure it is a good change. It’s really easy to tell if there are no side effects. It’s much harder when there are.
@rok-cesnovar: thank you for laying out the options! Your opinion matters. What is decided here will be our collective burden once it’s in the code base.
@wds15: did Rok characterize the options correctly? Do you have a different preference?

I completely agree with you that we need to automate this. I thought I asked you for help in designing this, but if they wasn’t clear, please help. We should do this sooner rather than later.

I think this looks really promising and that at least makes our end-to-end performance tests semi-automatic. Which is probably enough as long as we add additional test models to have better performance test coverage

For the updates to the Boost library, do we need to do more here? Don’t we need to tell people we have changed our implementation to match a higher boost version? Should this be in the upgrade boost script just in case we upgrade boost to a version without this patch?

I dont think we need to do more. For 2 reasons:

Sebastian has added a test here that will fail if Boost with no patch is used
The next Boost version that will actually be eligible for upgrade (the next Boost version that will have Boost Headers package on CRAN) will have the patch already done

I dont think we need to do more.

That's a good enough justification for me!

@wds15: we're waiting on a few comment clarifications. I tried rereading your comments (in the code) and I don't understand what's being written which means I won't understand what's there a few months from now.
Once that's cleaned up, I'll approve and we can merge.

Great. I will handle those on the weekend. Thanks.

Review comments should be adressed now. The comments are now more explanatory and are hopefully  clear by now.
I don't think we need a new issue.
(though I will try out once if we can avoid the double to long double propagation of boost::math functions for all of stan::math by using an adequate define)

(though I will try out once if we can avoid the double to long double propagation of boost::math functions for all of stan::math by using an adequate define)

This should be a separate issue.  Is there a separate way to configure this with some kind of environment variable for the build rather than hard coding it?  I wouldn't want to make the decision dynamic as that would introduce more overhead.

To be clear, the global precision thing for boost will be a different issue, of course.

@wds15 I resolved a simple conflict that came after the merge of #1340 . Hope that is fine with you.

Thanks @rok-cesnovar .... now there is another hickup which looks new to me. Probably to some more tighter style things? If you see whats going on then any hand is great.

Oh, I missed that, will fix, sorry.

It was just a bad job on my end. Fixed now.

BTW - requiring that issues are created before approving does convey a serious portion of mistrust. At least this is the impression I am getting out of this. I don't think that this is a good way of collaborating.

Can you elaborate more about this? Happy to take it offline if you'd like to speak more frankly.
It's not meant to be mistrust. It's pragmatic using the tools we have on hand, i.e. GitHub.
Maybe we should rotate who does release notes.... using GitHub, until recently, it was very difficult to get a list of PRs (they've changed that). So to handle it, we went with issues describing what's happening in our code base and PRs handling the actions. Having the issues there before merging PRs allow us to close issues with the PRs.

A mere, "please create it" would be sufficient.

I think the timing is the issue? You'd rather I requested you create an issue and merge this PR first?
That seems like bad practice. It's a requirement of getting the PR in: it's part of the documentation of this particular issue. As a policy for everyone, issues must be created prior to the PR getting in. Let's just not get in the habit of making exceptions... the more established devs should be following these guidelines and setting examples rather than looking to bend the guidelines.

On a separate note, I'm editing the top-level comment to include the new issue so we can easily know this PR closed 2 issues.

I've always been a little confused about the issue thing. What does an
issue provide that a PR doesn't? In GitHub, PRs are issues. Not saying
the policy is bad but just wondering what the goals are and how to think
about it.
In my past life as a professional engineer an issue would describe a user's
perspective and a PR a technical one. So for technical issues that don't
affect users we wouldn't have a separate issue.
Here it seems like an issue can say the exact same thing as a PR.
For release notes, our current system is really bad as neither the issue
nor the PR tend to be written from a user's perspective, especially not a
Stan user (vs Math user). One thing I like at Google is that every PR has a
section for what to include in the release notes for if that PR is merged.
On Sat, Sep 14, 2019 at 08:52 Daniel Lee<EMAIL_ADDRESS>wrote:

BTW - requiring that issues are created before approving does convey a
serious portion of mistrust. At least this is the impression I am getting
out of this. I don't think that this is a good way of collaborating.
Can you elaborate more about this? Happy to take it offline if you'd like
to speak more frankly.
It's not meant to be mistrust. It's pragmatic using the tools we have on
hand, i.e. GitHub.
Maybe we should rotate who does release notes.... using GitHub, until
recently, it was very difficult to get a list of PRs (they've changed
that). So to handle it, we went with issues describing what's happening in
our code base and PRs handling the actions. Having the issues there before
merging PRs allow us to close issues with the PRs.
A mere, "please create it" would be sufficient.
I think the timing is the issue? You'd rather I requested you create an
issue and merge this PR first?
That seems like bad practice. It's a requirement of getting the PR in:
it's part of the documentation of this particular issue. As a policy for
everyone, issues must be created prior to the PR getting in. Let's just not
get in the habit of making exceptions... the more established devs should
be following these guidelines and setting examples rather than looking to
bend the guidelines.
On a separate note, I'm editing the top-level comment to include the new
issue so we can easily know this PR closed 2 issues.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAGET3BT2GIPNIPMEPNWRRTQJTNCJA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD6W3IDA#issuecomment-531477516,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAGET3ELKFYHTNCQOPNXSLTQJTNCJANCNFSM4HRGN7BA
.

I've always been a little confused about the issue thing. What does an issue provide that a PR doesn't? In GitHub, PRs are issues. Not saying the policy is bad but just wondering what the goals are and how to think about it.

A technical note: in GitHub, PRs are not issues. They started out as separate entities, but now they're a lot closer together with similar functionality. (There's some history here, Why is it that PUll Requests and Issues share the same numbering, but it's not really that useful.)
Present day, they share a lot of the same fields, but there are a few notable differences:

issues don't have code. PRs do.
PRs can close issues. Issues can't close PRs. See: closing issues using keywords

In my past life as a professional engineer an issue would describe a user's perspective and a PR a technical one. So for technical issues that don't affect users we wouldn't have a separate issue.

We're not too far off from what you're talking about. Ideally, the issue describes "what" and the PR describes "how." The users of the Math library are C++ developers, so knowing what sorts of changes are happening in the codebase is sufficient. In order to get a sense of what details are involved is the how. For example, an issue would state something like "stan/math/fwd/mat/meta/operands_and_partials.hpp has incorrect includes" and the body of the issue would describe what's wrong. The PR would address the code to fix that issue (or multiple issues) and it would have all the technical detail that's associated with the fix. If that requires some design discussion, that's in the PR, not really in the issue unless that's the point of the issue.

For release notes, our current system is really bad as neither the issue nor the PR tend to be written from a user's perspective, especially not a Stan user (vs Math user). One thing I like at Google is that every PR has a section for what to include in the release notes for if that PR is merged.

I agree: it's not great. Our issues need a lot of work curating. Ideally, our release notes process should just be... take the titles from the issues that were closed by PRs since the last release. That should tell the story of how the codebase is changing.
For what you're requesting, something that's designed for a Stan user, that's something different and we won't have an automated process for that. If you think of something, please post it.

On Tue, Sep 24, 2019 at 01:32 Daniel Lee<EMAIL_ADDRESS>wrote:

I've always been a little confused about the issue thing. What does an
issue provide that a PR doesn't? In GitHub, PRs are issues. Not saying
the policy is bad but just wondering what the goals are and how to think
about it.
A technical note: in GitHub, PRs are not issues. They started out as
separate entities, but now they're a lot closer together with similar
functionality. (There's some history here, Why is it that PUll Requests
and Issues share the same numbering
https://github.community/t5/How-to-use-Git-and-GitHub/Why-is-it-that-Pull-Requests-and-Issues-share-the-same-numbering/td-p/22928,
but it's not really that useful.)
Present day, they share a lot of the same fields, but there are a few
notable differences:

issues don't have code. PRs do.
PRs can close issues. Issues can't close PRs. See: closing issues
using keywords
https://help.github.com/en/articles/closing-issues-using-keywords

That first citation is a doozy :) anyway, you're right that a PR has code
attached. My point is that they're numbered the same way and PRs have all
the same explanatory text fields as an issue.

https://help.github.com/en/articles/closing-issues-using-keywords

In my past life as a professional engineer an issue would describe a
user's perspective and a PR a technical one. So for technical issues that
don't affect users we wouldn't have a separate issue.
We're not too far off from what you're talking about. Ideally, the issue
describes "what" and the PR describes "how." The users of the Math library
are C++ developers, so knowing what sorts of changes are happening in the
codebase is sufficient.

I think this is the sentence where this completely diverges from what I was
describing (not that we have to aim for what I was describing). I would
first of all note that the vast overwhelming majority of the users of the
Math library are not developers but are Stan users.
But for the sake of describing what I'm talking about a little further,
let's aim for math library developer users as you say. Describing code
changes in the issue text is not what my previous employers would have
wanted in this dichotomy. Math users want to know how the issue affects
their use of the library. If include file ordering causes a bug or error
message, the issue would describe that user-facing symptom. If refactoring
serves some use case, the issue would describe that. Even the purely
internal technical issues we have that Math users would not perceive (e.g.
the flattening) could be described in terms of its high level goals instead
of the technical steps required to accomplish those implicit goals.
This may best be illustrated by example.

Adding a new solver (or any function) is easy - issue described how to
use the solves and what it's useful for.
If you know ahead of time that some refactoring is needed for performance
improvements, you can describe the characteristic of those improvements and
why the technical change will achieve those. Other performance improvements
could be discovered serendipitously and have no associated issue.
A refactor to allow multiple translation units would have an issue
describing how and why folks use multiple translation units before possibly
talking about the multiple strategies you could use to achieve that.

In order to get a sense of what details are involved is the how. For

example, an issue would state something like "stan/math/fwd/mat/meta/operands_and_partials.hpp
has incorrect includes" <#m_-9137811949223574004_1253> and the body of
the issue would describe what's wrong. The PR would address the code to fix
that issue (or multiple issues) and it would have all the technical detail
that's associated with the fix. If that requires some design discussion,
that's in the PR, not really in the issue unless that's the point of the
issue.
For release notes, our current system is really bad as neither the issue
nor the PR tend to be written from a user's perspective, especially not a
Stan user (vs Math user). One thing I like at Google is that every PR has a
section for what to include in the release notes for if that PR is merged.
I agree: it's not great. Our issues need a lot of work curating. Ideally,
our release notes process should just be... take the titles from the issues
that were closed by PRs since the last release. That should tell the story
of how the codebase is changing.

Well as you note we have two sets of users. If Math is focused on writing
issues and PRs for an audience of Math users, this will obviously never
work.
Would you be comfortable switching that prescribed audience for Math? An
alternative would be for conscientious devs to include a release notes
section in the issue.
For what you're requesting, something that's designed for a Stan user,

that's something different and we won't have an automated process for that.
If you think of something, please post it.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/stan-dev/math/pull/1255?email_source=notifications&email_token=AAGET3DDWRQ75GGJYN7ATBLQLGQ6VA5CNFSM4HRGN7BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7NDX2A#issuecomment-534395880,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAGET3CKNEBQLYRPXTGN7F3QLGQ6VANCNFSM4HRGN7BA
.

But for the sake of describing what I'm talking about a little further,
let's aim for math library developer users as you say. Describing code
changes in the issue text is not what my previous employers would have
wanted in this dichotomy. Math users want to know how the issue affects
their use of the library.

By user of the math lib, I'm taking you to mean anyone who uses Stan, because it's all built on top of the math lib.
Do you have examples of release notes from these previous employers we could check out?  I'd like to see what you consider good examples.  Here are TFP's (you have to go past the first one, which has only very high-level notes).
https://github.com/tensorflow/probability/releases
I'm not a TFP user or dev, but it looks like a mix of user-facing and dev-facing stuff.
I don't think anyone's suggesting code changes in the issue text.   Daniel said:

Ideally, the issue describes "what" and the PR describes "how."

I think the trick's going to get the dev-facing issue/PR notes in Stan user-comprehensible form.  If we speed up binomial_logit_glm, a user won't necessarily know where that gets used in something like brms.
You can see from our current release notes that we're not even close to being end-user facing:

https://github.com/stan-dev/math/releases/tag/v2.20.0

The release notes in the announcements category of Discourse isn't comprehensive, though it does pull out some notes that are of interest to users:

https://discourse.mc-stan.org/t/stan-2-20-0-released/9797

It'd be great if we could describe the effects of changes in language that end users of Stan could understand.   For instance, I could've included those notes on compile times in the PR for that.
The other issue is how to get release notes in front of users.  Are  they available anywhere other than on the GitHub releases?

No. I was thinking of starting a tradition of the TWGD writing a release
notes summary for each release with the most salient Stan user-facing
changes and posting this on Discourse (and maybe also the Stan Github).

Daniel and I used to write these.  It'd be great if someone could take
the time to translate for the users.  I don't know how comprehensive they
should be, which is why I'm hoping to see some good open source examples.

Bob

If you wanted to some of our previous release notes: https://github.com/stan-dev/stan/releases?after=v1.3.0
Automation is good, but so is taking a step back and summarizing what's happened since the last release. We should try to do a mix that's reasonable.
And... yes, if you have other examples, please post them.

Why are you guys looking for examples of open source release notes? I'm confused. I thought this conversation was about issues vs. PRs.

Either way, @syclik it seems like you missed my proposal to you in the past few exchanges. The proposal: ask that Math issues be written from the perspective of the majority of users of the Math library (currently, that would be users from the Stan modeling language). What do you think about that?

No. I was thinking of starting a tradition of the TWGD writing a release
notes summary for each release with the most salient Stan user-facing
changes and posting this on Discourse (and maybe also the Stan Github).

Daniel and I used to write these.  It'd be great if someone could take
the time to translate for the users.  I don't know how comprehensive they
should be, which is why I'm hoping to see some good open source examples.

I think the PR release notes would make that a bit easier
https://github.com/stan-dev/math/issues/1357

On Sep 25, 2019, at 2:23 PM, seantalts<EMAIL_ADDRESS>wrote:
Why are you guys looking for examples of open source release notes? I'm confused. I thought this conversation was about issues vs. PRs.

You brought it up when you wrote:

One thing I like at Google is that every PR has a section for what to include in the release notes for if that PR is merged.

Daniel replied:

Ideally, our release notes process should just be... take the titles from the issues that were closed by PRs since the last release.

This entangled the issue of PRs vs. issues with how release notes are generated.
I'm confused about how this discussion arose on a merged PR.  I stopped trying to regulate where discussions happen when people started telling me it was counterproductive to be so process oriented.  Certainly less work for me just responding on-channel to whatever comes up.

@seantalts

The proposal: ask that Math issues be written from the perspective of the majority of users of the Math library (currently, that would be users from the Stan modeling language). What do you think about that?

I don't think what you have written works. A modification might be that each issue includes something like that.
One role issues have is to break down work into pieces that can be worked on. Not all of that fits nicely into a form that's useful for Stan users.
I think what @rok-cesnovar has proposed is good: take the issues written in Math for the Math audience (a mix of direct Math users, Stan users, and mostly Math devs) and then write Stan-user-facing notes after this is done.

This is what he proposed on Discourse: Stan Math Release notes draft