Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Layer pybindings #20

Open
wants to merge 25 commits into
base: master
Choose a base branch
from
Open

Conversation

atrettin
Copy link
Contributor

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 atrettin marked this pull request as draft August 31, 2020 14:24
@atrettin atrettin marked this pull request as ready for review September 17, 2020 14:42
@atrettin
Copy link
Contributor Author

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
Copy link
Contributor Author

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.

This commit uses the low-pass filtering introduced in the latest PR to SQuIDS. This low-pass filtering can be applied while integrating the interaction picture states over time by calling
`Set_EvolLowPassCutoff` to choose a cut-off frequency of the filter. The frequency has units 1/length. This filter can speed up calculation of interaction picture states significantly in the presence of
eV-scale sterile neutrinos.

There are a few more additions to the Python interface of `nuSQUIDSLayers` that allow the filter to be applied also when projecting out probabilities from interpolated states.
… mode

Pybindings are greatly improved with overloads instead of differently named
functions, with doc-strings for the most important methods.

Preparations have been made to enable `both` mode to propagate neutrinos
and antineutrinos in the future. Work to enable it properly are still
ongoing.
The single-energy mode of nuSQuIDS is now initialized in the same way
as for the multiple-energy mode. This enables propagation of neutrinos
and antineutrinos at the same time.

Pybindings for nusquids layers have been improved further to allow
passing of different averaging ranges for every event.
@arguelles arguelles self-requested a review December 24, 2020 21:46
@arguelles arguelles self-assigned this Dec 24, 2020
@arguelles
Copy link
Owner

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
Copy link

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
Copy link
Collaborator

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
Copy link

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants