adding an exponential disk + Hernquist bulge model

[CarrieFilion]

Added the option to make an exponential disk + Hernquist bulge cylindrical basis. I did a few minor sanity checks and would be happy to run any additional checks that you can think of! Also happy to change naming etc if there is a specific convention that is desired.

Note - in the Hernquist portion, I transform the cylindrical R and z into spherical rr for the density profile. Not sure if this is the right thing to do

I also found some minor typo issues in BiorthBasis.cc for the double exponential model.

[The9Cat]

Just to make sure I understand: you want a cylindrical basis that has a biorthogonal basis whose lowest order basis function that is the sum of a exponential basis and Hernquist-model bulge? So the idea is that you want to simulate a single component disk with this new profile?

[The9Cat]

Is there some reason that this can not be accommodated by adding a separate bulge component?

[CarrieFilion]

Indeed - I would like the lowest order basis function to be the sum of a disk and a bulge. I do not intend to run simulations with this basis, I want to use this basis to analyze cosmological simulations. I had initially started off with separate disks and bulges, but then I would need to assign stars to the disk or to the bulge which created a new set of problems. I figured that a single basis would be more straightforward - I can just determine the disk + bulge properties from density fits and describe the total stellar potential with this one basis (Mike also suggested this!)

[michael-petersen]

Nice changes! I think just a little tweak to make.

[michael-petersen]

I think this line should be

Suggested change

	cxxopts::value<string>(dmodel)->default_value("0.1"))
	cxxopts::value<string>(HERNA)->default_value("0.1"))

[michael-petersen]

And this line should be

Suggested change

	cxxopts::value<string>(dmodel)->default_value("0.1"))
	cxxopts::value<string>(Mfac)->default_value("1.0"))

(both the variable name and the default value)

[michael-petersen]

We do offer the construction of custom-density-function bases, so I think it makes sense for us to also add some clearer documentation somewhere of how a user would do this, particularly since the changes can/should happen in a few places if one wants the density model fully available (as the edits here demonstrate!).

[The9Cat]

Yes, I would suggest that rather than change the code for this customized application that we provide a user-define density/potential functor for pyEXP analysis.

[The9Cat]

Sigh, looks like the DiskDens functor didn't make the cut on the last feature branch. I'd suggest that we restore that. It would not be much work; one simply needs to register all call out function for Cylindrical in BiorthBasis. Mike, can you recall why we axed that?

[UPDATE: I think I axed it because I had too much on my plate to gin up and example and document it. It would be straightforward to restore it. Maybe Carrie can try putting in her new joint function as test case?]

[michael-petersen]

As I see it, the main reason to keep this functionality (e.g. edits in cylcache and DiskModels.H) is to enable more powerful basis creation -- that is, not needing mpi4py to construct high-resolution (or many) bases.

That said, a functor in pyEXP would be a very flexible option for experimentation, and for the more patient. So I agree with that strategy in the long term.

[The9Cat]

As I see it, the main reason to keep this functionality (e.g. edits in cylcache and DiskModels.H) is to enable more powerful basis creation -- that is, not needing mpi4py to construct high-resolution (or many) bases.

True. But it doesn't solve the problem of how a user adds a new density target 'on the fly'. The proposed one is general useful, but suppose some particular project needs a more specialized target; is the only option to add it like this one, with a lists of new parameters, etc.? Is there a way to do this in a more modular way? As I scan the diffs, I note that most of this is about adding a lot of new parameters. The main problem with EXP, imho, is that we already have too many parameters. My bad but still...

An alternative to your mpi4py avoidance is to embed the Python interpreter with pybind11 and let the user supply the Python code for the new target density as a file or string. This does make provenance harder (unless you cache the Python density function itself, which would be easy enough). And the calls to the density function will be slower than native C++ but target is only called when computing the EOF covariance matrix. That's not the computational bottleneck, I don't think.

Anyway, feel free to proceed with this. I'm not trying to impede progress, but I am trying to learn from this PR if there is something we could be or should be doing better. It does bring up some interesting questions for future user needs. I will add the functionoid call from Python (it's pretty easy) after we merge #101.

Assuming that we are going ahead with this, can I ask if @CarrieFilion has checked whether you get a high quality basis out of this? Are there particular parameter choices that work or don't work?

[The9Cat]

Okay, here's another related idea: write a function class in exputils that binds a Python function to C++. Expose this through pybind11 to pyEXP and include a native C++ constructor as well. So the user provides a Python module with the desired disk_density function. Then the problem is solved on both sides once and for all. I wrote a simple example. Seems to work just fine.

[The9Cat]

I have set of changes to this PR that implements Python functors in all the same places as the diskbulge. I have not yet tested it (it does compile). This could be easily added (and probably should be added if we decide to keep this) to Cylinder.

[UPDATE: Functor implemented in Cylinder.]

[The9Cat]

Summary of the last commit blob

Fix a few minor issues identified while testing the user-specified Python density target

1. Remove duplicated diskbulge parameters from cylcache option parsing
2. Use an existing Python interpreting in global namespace if it's already created (e.g. by Jupyter)
3. Documentation for the form of the Python density function; also see the example in github/Minor updates pyEXP-examples#13 in the Tutorials/Basis directory.

TL;DR; usage

1. Use python as the dtype in your config or as a parameter for cylcache.
2. Supply the name of a Python module, pyname that implements a function called disk_density(R, z, phi).

For example, put the code below in the file pyDens.py and pass the parameter pyname: pyDens

import math

# Gaussian blob density
def disk_density(R, z, phi):
    a = 0.01                      # Scale radius
    h = 0.002                     # Scale height
    f = R*math.exp(-R*R/(a*a))
    g = math.exp(-z*z/(h*h))
    return f * g;

Comments

1. The code design allows for user specification of the function name (fixed by default to disk_density but I have not exposed that.
2. The EmpCylSL cache does not know anything about the intent of the target density, only the grid features. For example, if one made two bases with the same grid parameters but different target densities, the caches could be used interchangeably. Could be confusing. One possible fix is to name the target density function (e.g. the associated Dtype string) and EmpCylSL can require that the names match. However, Cylinder in exp relies on the user to specify a cache in common usage and only has the exponential target (and now Python) for standalone use. So maybe I'm not inclined to further mess with cache requirements. The user will have to be careful...

Tests

• Verify that cylcache makes a basis using the python method
• Verify that the pyEXP.basis.Basis.factory makes and reads a EmpCylSL cache
• Verity that Cylinder in exp makes and reads a EmpCylSL cache
• Verify that gendisk makes a basis and ICs
• Basis creation using Python functor works with MPI?
• Check diskbulge basis for sanity and orthogonality for astronomically realistic parameters

[michael-petersen]

This is definitely a good update! I have been working with the new dtype: python and pyname functionality, bases are turning out as expected (i.e. a exponential + sech^2 disc is the same basis as the native version).

I totally agree that the proliferation of options, even in a siloed piece of software like cylcache is not ideal. But now I'm wondering if we even want cylcache as extra maintainence, or if we should go 'all-in' on Python creation?

The basis creation is obviously quite computationally expensive on a single processor; I think we'll want big flashing warnings to users if they run an example notebook. With mpi4py the wallclock time is manageable, so pumping up that option is probably a good strategy. I also agree that the main shortcoming right now with the added flexibility is the possibility of losing the function definition for basis creation. I'm not sure what to do about this other than warn users to look carefully (printing the function is an option, but only if the user is going to look at it!).

All changes compile and work for me, so I'll approve this; we can think about deprecating cylcache in future versions. As for workable diskbulge parameters, I think users can always shoot themselves in the foot with non-ideal parameters, so is it a hard requirement for merging?

[The9Cat]

Thanks for checking all of this, Mike! This trick might come in useful in other cases.

I agree with you that cylcache seems like duplication. One can run basis creation from Jupyter using ipyparallel and mpi4py. Cell magic %%px is useful for that. Maybe it's a matter of taste? Along the same lines, cylcache could become a standalone Python app. Still, it's another executable but there would be 100% consistency in the interface. And most users would feel comfortable tweaking that. I'm in favor of deprecation, but not sure yet which direction forward is best.

I'm not overly concerned about sane parameters for diskbulge, but if the diskbulge scale (HERNA) gets large compared to the scale length, there could be grid mapping issues. Also, there could be basis conditioning issues. E.g. suppose one conditions with the exponential deprojection but has a massive small HERNA value relative to the disk scale length. Could get weird. But this will be the case with any user-specified target and parameter set, as you say.

[The9Cat]

One more comment: we have a bit of a cache file ambiguity even without the python target density. Stuff like gendisk and cylcache can use a variety of targets but EmpCylSL doesn't know about the function it's been passed. exp's Cylinder class has only one fallback target (exponential) but, of course, it now has python.

[The9Cat]

Let's wait for @CarrieFilion before merging. Also, there is the issue of which order to merge this one and #101; there will probably be conflicts, albeit easy to fix.