Lambda enhancements

[themikelau]

This pull request is to fix the issue of discontinuous Nanjing lambda's (Issue #418).

• Where the the radius R exceeds the maximum radius R_max for which the Nanjing polynomial fits work well (this maximum radius was ported over from Startrack), simply use the lambda at R_max, i.e. lambda(R>R_max) = lambda(R=R_max).
• Nanjing lambda's are now linearly interpolated between the masses (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 50, 100 Msun) and linearly in logged metallicities (Z = 0.001 'pop. II', Z = 0.02 'pop. I'). And where the input mass or metallicity falls out of its domain, we do not extrapolate linearly but merely use calculate lambda with the values at the domain boundary. E.g., lambda(Z, M > 100 Msun) = lambda(Z, M = 100 Msun); lambda(Z < 0.001, M) = lambda(Z < 0.001, M).
• Lambda's now take as an input the effective ZAMS mass, m_Mass0, (ZAMS mass corresponding to rejuvenated star?) rather than the true birth mass. Note that this is opposite to Startrack's implementation.
• The purely gravitational lambda is now capped to be below one, lambdaG <= 1. lambdaG exceeding one is unphysical, as lambdaG = 1 corresponds to the extreme limit of all envelope mass being located at the stellar radius.
• The lambda with internal energy is now constrained to exceed the purely gravitational lambda, lambdaB >= lambdaG. This is because internal energy must act to make the envelope easier to eject.

In terms of code, here's a summary of what has been done to achieve the above.

• Change the CalculateLambdaNanjing functions in CHeB.cpp, EAGB.cpp, TPAGB.cpp, and HG.cpp to take mass and metallicity as inputs to facilitate interpolation, e.g. double CHeB::CalculateLambdaNanjing(double mass, double metallicity) {...}.
• Add a double BaseStar::CalculateMassInterpolatedLambdaNanjing(const double metallicity) that performs mass interpolation for given metallicity. This function accesses the star's m_Mass0, calls CalculateLambdaNanjing for its adjacent mass bins, then linearly interpolate between the two values.
• Add a double BaseStar::CalculateMassAndZInterpolatedLambdaNanjing() that performs metallicity interpolation. This function accesses the star's metallicity, calls CalculateMassInterpolatedLambdaNanjing for its adjacent metallicity bins, then linearly interpolate between the two logged metallicities.
• When evaluating lambda for a given radius, use the value min(Rmax, m_Radius) instead of m_Radius, where Rmax is different for each mass.
• Capped lambda by adding in CalculateLambdaNanjing:

lambdaBG[1] = std::min( std::max(0.05, lambdaBG[1]), std::min(1.0, maxBG[1]) );         // clamp lambda G to [0.05, min(1,maxG)]
lambdaBG[0] = std::max( std::min(lambdaBG[0],maxBG[0]), std::max(0.05, lambdaBG[1]) );  // clamp lambda B to [ max(0.05,lambdaG), maxB]

What could be improved

• This needs to be coded up in a way that would allow a switch between the enhanced and old lambdas to be easily implemented.
• double Rmax = 1000000.0; is ugly?
• Would it be better to have CalculateLambdaNanjing take in discrete mass and metallicity values, since CalculateMassAndZInterpolatedLambdaNanjing and CalculateMassInterpolatedLambdaNanjing are doing all the interpolation anyway? E.g. Having something like switch (mass) { case '1.0': [...] case '2.0': [...] instead of if mass < 1.5 [...] elseif mass < 2.5 [...]?

Overall, this is all very rudimentary I would appreciate feedback from code experts. Once the code is polished up, I will update this branch with the latest updates from dev (there are currently conflicts) and update the changelog.

[themikelau]

I need to fix a bug related to the mass and metallicity bins... Should be using M = [1.0, 2.0, ...] instead of M = [1.5, 2.5, ...]

[jeffriley]

Overall, this is all very rudimentary I would appreciate feedback from code experts. Once the code is polished up,
I will update this branch with the latest updates from dev (there are currently conflicts) and update the changelog.

Good job. Just one thing: where you have parameters to functions (e.g. double TPAGB::CalculateLambdaNanjing(double mass, double metallicity) and others), our convention is (see COMPAS doc, ~p53):

(a) class member variable names begin with m_, and are typically followed by an uppercase alpha character (though that (the uppercase character) is sometimes relaxed if it makes sense - e.g. m_dM rather than m_DM)
(b) function parameter names begin p_ and, as for class member variables, are typically followed by an uppercase alpha character (with the same caveat)
(c) local variable names begin with a lowercase alpha character (no prefix as for class member variables and function parameters). Again, relaxed if it makes sense (e.g. maybe use ECSN rather than eCSN)
(d) function names should begin with an uppercase alpha character
(e) all names (function names, variable names) should be camelCase (rather than snake_case)

The rationale for (a)-(e) is consistency; the rationale for (a), (b), and (c) in particular is that readers of the code at some later stage can tell easily where a variable came from - whether it was calculated locally, passed in as a parameter, or whether it is a class member variable.

Typically we don't modify parameters that were passed in to functions - although there are a (very) few places in the code where we do for pragmatic reasons (performance mostly). If a function modifies either class member variables or parameters it should be documented in the descriptive text before the function declaration.

If parameters are not modified by the function they should be declared const, and similarly if the function does not modify class member variables it should be declared const.

double Rmax = 1000000.0; is ugly?

Yes, but maybe necessary. Could you use std::numeric_limits::max() instead of the arbitrary 1000000.0?

Would it be better to have CalculateLambdaNanjing take in discrete mass and metallicity values, since
CalculateMassAndZInterpolatedLambdaNanjing and CalculateMassInterpolatedLambdaNanjing are doing all the
interpolation anyway?
E.g. Having something like
switch (mass) { case '1.0': [...] case '2.0': [...] instead of if mass < 1.5 [...] elseif mass < 2.5 [...]?

I've only looked at the changes to the code through github - I haven't looked at the entire code in context -
so I'm not sure whether passing discrete values to CalculateLambdaNanjing() is better (I'll look later), but
you can't switch on a floating point number (and it would be a bad idea anyway - for the same reason we have
utils::Compare())

This needs to be coded up in a way that would allow a switch between the enhanced and old lambdas to be easily
implemented.

So, are we planning to make both the old and new lambdas available into the future? Or do we just want them available now for testing new vs old? Do we want to release a version on the code that allows users to choose, or do we just want to leave the old code in for a while so we can test?

If we want to release the code and allow users to choose between old and new lambdas, then we'll need a new program option. If we just want to leave both there for a while as we test, maybe just have a switch in the code that requires a code change - no
point coding up a new option if it's not permanent.

[themikelau]

Hi @jeffriley. Thanks for the detailed feedback! The requested changes related to variable names and type seem to be straightforward and I will fix them on Thursday.

so I'm not sure whether passing discrete values to CalculateLambdaNanjing() is better (I'll look later)

My reasoning is that in my changes, CalculateLambdaNanjing(...) no longer directly outputs a star's lambda, but has switched roles to become a lookup table for Xu & Li (2010) fitting parameters. And of course, this table has discrete mass and metallicity values since Xu & Li (2010) ran a discrete grid of stellar models. The new functions I have added, CalculateMassInterpolatedLambdaNanjing() and CalculateMassInterpolatedLambdaNanjing(const double metallicity), output lambda for continuous mass/metallicity by calling CalculateLambdaNanjing() at the discrete masses/metallicities and interpolating over them. So I thought that it would be odd and perhaps unsafe to allow CalculateLambdaNanjing() to take continuous mass or metallicity values.

So, are we planning to make both the old and new lambdas available into the future? Or do we just want them available now for testing new vs old? Do we want to release a version on the code that allows users to choose, or do we just want to leave the old code in for a while so we can test?

I think Ilya's suggestion has been to make a new program option for users to choose between old and new lambdas. I guess the purpose of this is so that the results of the old and new lambdas can be compared even with future versions of COMPAS.

[themikelau]

I have checked that the default code behaviour (not switching on any of the newly introduced lambda options) reproduces the detailed system plot (Fig 15 of methods paper). I also find identical output from 1000 binaries evolved using this new version and evolved using the existing version (v2.26.02).

[ilyamandel]

@themikelau : thank you very much for this PR!

It seems you updated the default pythonSubmit.py, but I think we have now changed to using runSubmit.py -- can you please update that instead?

Did you update the documentation to include the new code options?

A couple of comments on the code, though I'll defer to @jeffriley :

• I find CalculateLambdaNanjing() and EvaluateLambdaNanjing() slightly confusing; maybe rename the old CalculateLambdaNanjing into CalculateLambdaNanjingStarTrack, and the new EvaluateLamdaNanjing into CalculateLambdaNanjing?
• You often use multiple return statements inside complicated nested if/else loops. I think the usual COMPAS convention adopted by Jeff -- which I agree makes it easier to read the code -- is to assign the result to a variable and only return at the end of the function.
• (Perhaps for the future?) Instead of treating mass and metallicity interpolation separately, I might have opted for a single utility to interpolate on a 2-d grid (or not, if beyond the grid or if the interpolation along a given dimension is off). This might make for fewer special cases / functions. E.g., if I want the value at point {x,y} that falls within the grid, so that x_1<x<x_2 and y_1<y<y_2, and interpolation is on along both axes, then f(x,y)= dx dy f(x_1,y_1) + dx (1-dy) f(x_1,y_2) + (1-dx) dy f(x_2,y_1) + (1-dx) (1-dy) f(x_2,y_2), where dx=(x-x_1)/(x_2-x_1) and dy = (y-y_1)/(y_2-y_1).

[jeffriley]

You often use multiple return statements inside complicated nested if/else loops. I think the usual COMPAS convention adopted by Jeff -- which I agree makes it easier to read the code -- is to assign the result to a variable and only return at the end of the function.

Indeed. It's ok to do a sanity check or two at the very start of a function and, based on those sanity checks, return immediately if we shouldn't even execute the function, but after that, it's much easier for a future reader of the code to read and understand what's happening if they can assume there is one entry point (some languages allow multiple entry points for functions...) and one exit point for the function. So, as @ilyamandel says, the convention we use in COMPAS is to assign the return value to a variable, and return that at the end - and use if/else to execute only those parts of the function that are required rather than return early.

I'll repeat my mantra here. We should always code with future readers/maintainers of the code in mind. The compiler really doesn't care what our code looks like - it's very adept at parsing code however we format it. Code is typically written once, and read many times - if by being a little more verbose, formatting the code so it's easy on the eyes (even adding parentheses that aren't strictly necessary just for clarity for the reader), sometimes even repeating a little bit of code instead of trying to be tricky or get everything done in the fewest lines of code possible, we can make it easier to read and maintain (especially for people working on the code long after we've moved on), then that's what we should do. Compilers these days are almost always optimising compilers, so even though you might think a little bit of verbosity, or less trickiness, might make the code run slower, it probably won't - the compiler's going to optimise your code anyway...

[themikelau]

Thanks @ilyamandel and @jeffriley for your wisdom.

• I have opted for returning a single variable instead of using multiple return statements in the functions that calculate lambda.
• I have also updated the lambda options in compasConfigDefault.yaml (P.S. usage of runSubmit.py doesn't seem to be in the documentation.)
• I have renamed the old CalculateLambdaNanjing to CalculateLambdaNanjingStarTrack. And the new EvaluateLamdaNanjing wrapper function to CalculateLambdaNanjing.
• Updated documentation.
• Let's leave the 2-d interpolation function for the future.

Just to be safe I will run another test to check that the output with default options is same as in dev after these changes.

[themikelau]

Ok, I have just checked that the output of 1000 binaries is identical with that produced with v2.26.03. Happy for this to be merged.

[ilyamandel]

Re: runSubmit documentation -- @avigna is working on that, see #728 . The rest of the changes look great, thank you!