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

Add 4pi beam convolution #338

Draft
wants to merge 54 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
54 commits
Select commit Hold shift + click to select a range
2a71b78
rebuild the branch
mreineck Oct 30, 2024
3ab3f93
more comments
mreineck Oct 30, 2024
b4d05ed
update __init__.py
mreineck Oct 30, 2024
8494b5a
skeleton of the interface added
paganol Oct 31, 2024
b671b68
adjust error messages
mreineck Nov 4, 2024
8920c56
apply fix suggested by Luca Pagano
mreineck Nov 7, 2024
b56c6ea
add MuellerConvolver
mreineck Nov 7, 2024
40eceda
add standard 4pi convolution case
mreineck Nov 7, 2024
660d234
remove old, commented-out code
mreineck Nov 7, 2024
729c3cb
adjust docstring
mreineck Nov 13, 2024
834d1f2
Add dataclass BeamConvolutionParameters
ziotom78 Dec 6, 2024
e38f8fd
Merge branch 'master' into beam_comvolution_2
paganol Dec 6, 2024
98522ca
some bugs cleaned, now the convolution works
paganol Dec 6, 2024
4ebaf5b
Add functions to synthetize the a_ℓm of a Gaussian beam
ziotom78 Dec 6, 2024
e1be9ae
Add functions to synthetize the a_ℓm of a Gaussian beam
ziotom78 Dec 6, 2024
af6de79
method of simulation class added
paganol Dec 6, 2024
1a45ea3
Merge branch 'beam_comvolution_2' of https://github.com/litebird/lite…
paganol Dec 6, 2024
cf27ac2
some fixes
paganol Dec 6, 2024
4e34adc
4pi convolution now works
paganol Dec 7, 2024
6fc683d
Update ducc0 dependency in poetry.lock
ziotom78 Dec 10, 2024
cb77c00
Add new interactive notebook showcasing Gaussian beam synthesis
ziotom78 Dec 10, 2024
944b647
Merge remote-tracking branch 'origin/master' into beam_comvolution_2
ziotom78 Jan 9, 2025
e73830f
Add new dataclass SphericalHarmonics and make MBS return it
ziotom78 Jan 9, 2025
931de1d
Support nstokes=3 in SphericalHarmonics, reformat and add tests
ziotom78 Jan 10, 2025
17a9ea4
Add example notebook for beam synthesis
ziotom78 Jan 10, 2025
4046b52
Improve SphericalHarmonics and rename one method
ziotom78 Jan 14, 2025
e33575b
Add tests for the SphericalHarmonics constructor
ziotom78 Jan 14, 2025
7e8e732
Start updating beam_convolution so that it uses SphericalHarmonics
ziotom78 Jan 14, 2025
7d60b11
beam_convolution interface adjusted, some fixes on mbs, new function …
paganol Jan 16, 2025
ed515c2
automatic resize added
paganol Jan 16, 2025
55cb40a
comments in mbs removed
paganol Jan 16, 2025
c0d09ec
sky alms now adjusted automatically in add_convolved_sky_to_one_detector
paganol Jan 16, 2025
4ed706a
nthreads interfaced
paganol Jan 16, 2025
c76ccf6
interfaces updated
paganol Jan 16, 2025
d9ba1f5
fix bug
paganol Jan 16, 2025
831ee28
sign of hwp_angle inverted in mueller_convolver
paganol Jan 20, 2025
7a407a9
some fixes in blm generation
paganol Jan 22, 2025
1cbdacf
changes the equation of convolution
Jan 24, 2025
bd99c0c
change the sign of HWP angle
Jan 24, 2025
2c1667b
merge master
paganol Feb 6, 2025
dd9d86b
beam_convolution compatible with new pointing API
paganol Feb 6, 2025
13a22d7
placeholder for documentation added
paganol Feb 7, 2025
0426770
change of gauss_beam_to_alm parametrization
paganol Feb 7, 2025
713850c
test fixed
paganol Feb 7, 2025
3d69754
some fixes in beam_convolution
paganol Feb 17, 2025
26c7dca
Add `strict_typing` field to `BeamConvolutionParameters`
ziotom78 Feb 17, 2025
a36ae6f
generation of beam alms
paganol Feb 17, 2025
b5f0272
Merge branch 'beam_comvolution_2' of github.com:litebird/litebird_sim…
ziotom78 Feb 17, 2025
82b70b2
Fix typo in comment
ziotom78 Feb 17, 2025
9b031e0
Propagate `strict_typing` and add tests
ziotom78 Feb 17, 2025
96aabe9
Merge branch 'beam_comvolution_2' of github.com:litebird/litebird_sim…
ziotom78 Feb 17, 2025
d2370f5
interface of scan_map generalized
paganol Feb 17, 2025
76505b5
scan_map improved
paganol Feb 18, 2025
92d2dc1
simple notebook showing convoltuion added
paganol Feb 18, 2025
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .github/workflows/tests.yml
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@ jobs:
for proc in 1 5 9 ; do
echo "Running MPI test ($MPI) with $proc processes"
PYTHONPATH=. mpiexec -n $proc python3 ./test/test_mpi.py
PYTHONPATH=. mpiexec -n $proc python3 ./test/test_detector_blocks.py
done

- name: Tests OpenMPI
Expand All @@ -95,4 +96,5 @@ jobs:
for proc in 1 2 ; do
echo "Running MPI test ($MPI) with $proc processes"
PYTHONPATH=. mpiexec -n $proc python3 ./test/test_mpi.py
PYTHONPATH=. mpiexec -n $proc python3 ./test/test_detector_blocks.py
done
2 changes: 1 addition & 1 deletion .pre-commit-config.yaml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
repos:
- repo: https://github.com/astral-sh/ruff-pre-commit
# Ruff version.
rev: v0.1.7
rev: v0.1.15
hooks:
# Run the linter.
- id: ruff
Expand Down
75 changes: 75 additions & 0 deletions Beam synthesis example.ipynb

Large diffs are not rendered by default.

15 changes: 15 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,19 @@
# HEAD

- **Breaking change**: Change to the pointing API [#358](https://github.com/litebird/litebird_sim/pull/358), in detail:

1. DetectorInfo has three new attributes: pol_angle_rad (polarization angle), pol_efficiency (polarization efficiency)and mueller_hwp (mueller matrix of the HWP).

2. get_pointings() return only the orientation ψ of the detector, the polarization angle is a separate variable stored in the `Observation` class. The same class also handles the mueller_hwp for each detector, and it has a new bool variable `has_hwp` that is set to true if an HWP object is passed to `prepare_pointings()`.

3. The mock vPTEP IMo has been updated accordingly.

4. The `HWP` class has a new field called mueller, that contains the mueller matrix of the HWP.

5. The function `scan_map()` now handles three possible algebras: (i) no HWP, (ii) ideal HWP, (iii) generic optical chain.

- Implementation of distributing detectors across the MPI processes by grouping them according to given attributes [#334](https://github.com/litebird/litebird_sim/pull/334)

- **Breaking change**: `PointingSys()` now requires `Observation` as an argument. And several functions to add pointing systematics are merged into `left_multiply_syst_quats()`.

- **Breaking change**: Redefinition (ωt instead of 2ωt) of the hwp angle returned by `get_pointings()`. Change in the pointing returned by `Observation.get_pointings()`, now behaving as it was before this [commit](https://github.com/litebird/litebird_sim/pull/319/commits/b3bc3bb2049c152cc183d6cfc68f4598f5b93ec0). Documentation updated accordingly. [#340](https://github.com/litebird/litebird_sim/pull/340)
Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ Clicking on a dot reveals the `DetectorInfo` for that detector in real time, hig
Additionally, if you agree during the conversation to generate a detector file,
a list of starred detectors will be saved into a text file at the designated location after you closed the plot.

![plot_fp_usage](https://private-user-images.githubusercontent.com/83496454/388083441-93876c66-8b61-40c0-b23f-79a6998e5e33.gif?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MzIzNjAzOTAsIm5iZiI6MTczMjM2MDA5MCwicGF0aCI6Ii84MzQ5NjQ1NC8zODgwODM0NDEtOTM4NzZjNjYtOGI2MS00MGMwLWIyM2YtNzlhNjk5OGU1ZTMzLmdpZj9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDExMjMlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQxMTIzVDExMDgxMFomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWJkMmM3NTYzMzc4MTU2NGU3ODE2YWMzZTFlZDVjZmEzODJhNWMwMjhhMDY4MzUxMWU3OGYzZmZiYWFiZDcwNDQmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0In0.zG2QVKA9VeL5j79jTjiPNQ9kH8uyxVSeCzd3GBj3QVk)
![FP Usage](./images/plot_fp_usage.gif)

## Roadmap

Expand Down
Binary file modified default_imo/schema.json.gz
Binary file not shown.
31 changes: 31 additions & 0 deletions docs/source/beam_convolution.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
.. _beamconvolution:

Convolve Alms with a Beam to fill a TOD
=======================================

Minimal documentation should contain:
- quick descrition of the math (reference to papers)
- conventions
- details on the interface and examples
- description of the SphericalHarmonics class and how to store alms
- example with high level interface


The framework provides :func:`.add_convolved_sky`, a routine which
takes sky and beam alms perform harmonic space convolution and fills
the detector timestreams. It implements both the algebra for...

Methods of the Simulation class
-------------------------------

The class :class:`.Simulation` provides the function
:func:`.Simulation.convolve_sky`, which takes sky alms and perform the
convolution

API reference
-------------

.. automodule:: litebird_sim.beam_convolution
:members:
:undoc-members:
:show-inheritance:
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/images/detector_groups_case1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/images/detector_groups_case2.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/images/detector_groups_case3.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/images/grid_communicator.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/images/time_block_communicator.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
41 changes: 39 additions & 2 deletions docs/source/map_scanning.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,39 @@ Scanning a map to fill a TOD

The framework provides :func:`.scan_map`, a routine which scans an
input map accordingly to the scanning strategy and fills the detector
timestreams. You can fill with signal an existing TOD by using the
timestreams. It implements three possible algebras:

- No HWP:

.. math::
d_t = T + \gamma \left(Q\cos(2\theta + 2\psi_t)+U\sin(2\theta + 2\psi_t)\right),
:label: noHWP

where :math:`\theta` is the polarization angle of the detecotor, :math:`\psi_t`
is the orientation of the detector at the time :math:`t`, and :math:`\gamma`
is the polarization efficiency.

- Ideal HWP:

.. math::
d_t = T + \gamma \left(Q\cos(4\alpha_t - 2\theta + 2\psi_t)+U\sin(4\alpha_t - 2\theta + 2\psi_t)\right),
:label: idealHWP

where :math:`\alpha_t` is the HWP angle at the time :math:`t`.

- Generic HWP:

.. math::
d_t = (1,0,0,0) \times M_{\rm pol} \times R(\theta) \times R^T(\alpha_t) \times M_{\rm HWP} \times R(\alpha_t) \times R(\psi_t) \times \vec{S}
:label: genericHWP

where
* :math:`M_{\rm pol}` is mueller matrix of the polarimeter;
* :math:`M_{\rm HWP}` is mueller matrix of the HWP;
* :math:`R` is rotation matrix;
* :math:`\vec{S}` is the Stokes vector.

You can fill with signal an existing TOD by using the
function :func:`.scan_map_in_observations`, as the following example
shows:

Expand Down Expand Up @@ -91,6 +123,7 @@ shows:
-0.00601
-0.00676

The code automatically selects the fastest algebra based on the provided HWP.

The input maps to scan can be either included in a dictionary with the name of
the channel or the name of the dectector as keyword (the routines described in
Expand All @@ -103,13 +136,17 @@ coherent, so either a single Observation and a single numpy array, or same
lenght list of Observations and numpy arrays.
If the input map is ecliptic coordinates set `input_map_in_galactic` to `False`.
The effect of a possible HWP is included in the pointing information, see
:ref:`scanning-strategy`.
:ref:`scanning-strategy`, the polarization angle of the detectors is taken from
the corresponding attributes included in the observations. The same applies to
the polarization efficiency.

The routine provides an on-the-fly interpolation of the input maps. This option
is available through the argument `interpolation` which specifies the type of TOD
interpolation ("" for no interpolation, "linear" for linear interpolation).
Default: no interpolation.

The low level function, :func:`.scan_map`, allows a more refined handling of the
inputs.

Methods of the Simulation class
-------------------------------
Expand Down
25 changes: 25 additions & 0 deletions docs/source/mpi.rst
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,31 @@ variable :data:`.MPI_ENABLED`::
To ensure that your code uses MPI in the proper way, you should always
use :data:`.MPI_COMM_WORLD` instead of importing ``mpi4py`` directly.

The simulation framework also provides a global object
:data:`.MPI_COMM_GRID`. It has two attributes:

- ``COMM_OBS_GRID``: This is an MPI communicator that contains all the
MPI processes with the global rank less than ``n_blocks_time * n_blocks_det``.
It provides a safety net to the operations and MPI communications
that are needed to be performed only on the partition of :data:`.MPI_COMM_WORLD`
that contain non-zero number of pointings and TODs. By default,
``COMM_OBS_GRID`` points to the global MPI communicator :data:`.MPI_COMM_WORLD`.
It is updated once :class:`.Observation` are defined. For example,
consider the case when a user runs the simulation with 10 MPI
processes but due some specific ``det_blocks_attributes`` argument
in :class:`.Observation` class, the number of detector and time
blocks are determined to be 2 and 4 respectively. Then the
simulation framework will store the pointings and TODs only on
:math:`2\times4=8` MPI processes and the last two ranks of :data:`.MPI_COMM_WORLD`
will be left unused. Once this happens, ``COMM_OBS_GRID`` on first 8
ranks (rank 0 to 7) will point to the local sub-communicator
containing the processes with global rank 0 to 7. On the unused
ranks, it will simply point to the NULL communicator.
- ``COMM_NULL``: If :data:`.MPI_ENABLED` is ``True``, this object
points to a NULL MPI communicator (``mpi4py.MPI.COMM_NULL``).
Otherwise it is set to ``None``. The user should compare
``COMM_OBS_GRID`` with ``COMM_NULL`` on every MPI process in order
to avoid running a piece of code on unused MPI processes.

Enabling/disabling MPI
----------------------
Expand Down
Loading
Loading