Layer pybindings

[atrettin]

For my upcoming light sterile analysis at low energies, there are a number of technical challenges for which I needed more low-level control over nuSQuIDS accessible from Python. The new class nuSQUIDSLayers is not terribly useful on its own as it is intended to be used as a backend for a Python implementation in PISA (https://github.com/IceCubeOpenSource/pisa). Since this code is going to be used not only for me but by anyone using nuSQuIDS in PISA, I'd like to make this a PR into the main branch of nuSQuIDS and would kindly ask the maintainers to review the implementation. I'm not very familiar with C++ and there are a number of points where things could probably be solved more efficiently and/or nicely by someone with more skill than me...

The requirements I'm attempting to fulfill are:

• reproducibility of other PISA tools: In PISA, we have a common module to calculate the distances and densities through the Earth's layers. The class nuSQUIDSLayers is initialized with those pre-computed layers and the results can be directly compared to those of other tools already implemented (prob3, GLoBES).

• externalized state interpolation: By default, nuSQuIDS uses linear interpolation on the interaction picture states. I wanted to experiment with other interpolation methods and the best way to do that is to simply provide Pybindings to pass an interpolated state. For that, I implemented EvalWithState, which takes in the components of the interpolated state as a numpy array

• averaging of fast oscillations: I will need the ability to smooth the extremely fast oscillations coming from light steriles at energies below 100 GeV, and that again needs to be accessible from Python. Ideally, I would also like the ability to use the approximation in the RHS of the differential equations.

[atrettin]

One thing that is just not yet where I want it to be is the multithreading. I had to construct it in such a way that all nodes can only work on the same layer in lock-step, meaning that all nodes have to wait for the slowest node at each layer. I couldn't get it to work in such a way that all nodes are truly independent and loop over all layers. It would be awesome if someone could help me with that! :-) See comment in L2087 of nuSQuIDS.h

[atrettin]

Also, it would be neat if the EvalWithState function could be turned into something that supports numpy broadcasting. That way, one could evaluate interpolated probabilities in 1D and 2D (even ND) without loops, and it might actually be faster! But again, that is simply above my capabilities, especially with Boost Python.

[arguelles]

I think these commits are good. Some functions on the headers do not have documentation and it would be nice to add that. Unless @cnweaver or @jsalvado are against I will merge them in ~1-2 weeks.

[austinschneider]

Hey @atrettin, I have a suggestion related to an architecture choice in your PR that may also help with the multi-threading.
Right now you are using the ConstantDensity body class and providing that to each individual nuSQUIDS object, swapping the body object for the next as needed. Another way to do this would be to implement a LayeredDensity body class that contains the information about all layers (currently stored within the nuSQUIDSLayers class) and provide that as the body to each individual nuSQUIDS object. I think this squares more nicely with the current architecture choices and has the added benefit that your current approach to the multi-threading needs almost no modification to achieve the "entire evolution through all layers (in each thread)" that you mention in the comments.

[cnweaver]

It's definitely possible to fix the multithreading in the current scheme (just make the task for each thread be to do one propagation through all of its layers, swapping the body as needed), but I have a patch in the works to enable multi-layered bodies directly. The main reason this can't be trivially done with a new subclass of the current base Body is that if you know there are sharp boundaries (as between constant density layers) you really want to start and stop the ODE solver at those points (which is one of the virtues of atrettin's implementation above), rather than forcing its adaptive stepping to laboriously discover where each boundary is by trial and error (or perhaps worse, not notice a boundary and produce imprecise results). Doing this for a single body turns out to be somewhat more difficult than it looks because of the very confusing way our time variable ends up being tracked between SQuIDS and nuSQuIDS.

[austinschneider]

Thanks for the followup @cnweaver. I agree, it's definitely possible to do it without a new body class. I hadn't considered the issue with the adaptive step size of the solver. It will be good to have the patch for multi-layered bodies so these things can be done without manually restarting the solver so to speak.

[atrettin]

I will soon try to split this PR up so that we can merge the more straight-forward additions. At least some functionality can be split up this way relatively easily, like the evolution low-pass filter. I've been optimizing my analysis and found that maybe the way I'm using the low-pass filter at the evaluation step is still not quite ideal and that I might want to add a linear ramp to the original SQuIDS averaging method instead.

[atrettin]

I've been thinking... basically, this class I implemented is just a container with many nusquids objects, and in principle this is totally separate from the layered body problem. The nusquids objects inside the container might as well use the Earth at a particular coszen or any other body for that matter. In my mind it would make most sense if @cnweaver has an implementation for layered bodies that I make this class more generic, such that one can pass any Body to the nodes inside it, that would make its existence more justified. What are your opinions?
(I'm working in the meantime to break up this big PR further).