ligo.skymap issueshttps://git.ligo.org/lscsoft/ligo.skymap/-/issues2024-03-18T22:52:17Zhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/48Document review signoff and approval workflow in Testing section of manual2024-03-18T22:52:17ZLeo P. SingerDocument review signoff and approval workflow in Testing section of manualhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/47Describe scope of unit tests, link to most important unit tests in Testing se...2024-03-18T13:09:19ZLeo P. SingerDescribe scope of unit tests, link to most important unit tests in Testing section of manualhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/45Add a few more historical events to the acceptance tests2024-02-20T17:51:52ZTito Dal CantonAdd a few more historical events to the acceptance testsSome proposals:
* A couple BBH events with precise skymaps: S200224ca and S200311bg (especially the latter, clearly detectable in three detectors)
* A very "ringy" skymap (HL-only event): S190707q
* A single-detector skymap (I see S19093...Some proposals:
* A couple BBH events with precise skymaps: S200224ca and S200311bg (especially the latter, clearly detectable in three detectors)
* A very "ringy" skymap (HL-only event): S190707q
* A single-detector skymap (I see S190930t on GraceDB, though I forget what that ended up like in GWTC)https://git.ligo.org/lscsoft/ligo.skymap/-/issues/44ligo-skymap-plot-volume does not respect transparent flag for the entire png2024-02-06T02:03:08ZGeoffrey Moligo-skymap-plot-volume does not respect transparent flag for the entire pngWhen creating volume renderings of sky maps using `ligo-skymap-plot-volume`, the `--transparent` flag is off by default. This is only respected for the subplot of the distance marginalization posterior, instead of the entire figure. This...When creating volume renderings of sky maps using `ligo-skymap-plot-volume`, the `--transparent` flag is off by default. This is only respected for the subplot of the distance marginalization posterior, instead of the entire figure. This makes it difficult (or impossible) to read the axes labels (see, e.g., https://git.ligo.org/emfollow/gwcelery/-/issues/752#note_927465).
Here's an example of two plots made of the bayestar skymap from https://gracedb.ligo.org/superevents/S240109a/view/.
Non-transparent plot made with
```plaintext
ligo-skymap-plot-volume bayestar.multiorder.fits -o plot_vol_no_transparent.png
```
![image.png](/uploads/8b4054cea21258549e86d30d42b46048/image.png){width=286 height=278}
Transparent plot made with
```plaintext
ligo-skymap-plot-volume bayestar.multiorder.fits -o plot_vol_transparent.png --transparent
```
![image.png](/uploads/985aa55d35535356fa456960d37b366f/image.png){width=295 height=289}
I think the intended behaviour should be to make the background of the entire figure (i.e., including all subplots) white if the transparent flag is not turned on, so that the labels are readable in all situations.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/43Catalog crossmatch example runs into an error with units2024-01-24T21:42:05ZGeoffrey MoCatalog crossmatch example runs into an error with unitsRunning the exact example code provided in https://lscsoft.docs.ligo.org/ligo.skymap/postprocess/crossmatch.html, I run into this error:
`UnitConversionError: 'Unit("Mpc2")' is not a scaled version of 'Unit("Mpc")'` and
`UnitConversion...Running the exact example code provided in https://lscsoft.docs.ligo.org/ligo.skymap/postprocess/crossmatch.html, I run into this error:
`UnitConversionError: 'Unit("Mpc2")' is not a scaled version of 'Unit("Mpc")'` and
`UnitConversionError: 'Mpc2' (area) and 'Mpc' (length) are not convertible` .
Here is the full traceback:
```python
10000it [00:00, 41431.73it/s]
---------------------------------------------------------------------------
UnitConversionError Traceback (most recent call last)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:982, in Quantity.to_value(self, unit, equivalencies)
981 try:
--> 982 scale = self.unit._to(unit)
983 except Exception:
984 # Short-cut failed; try default (maybe equivalencies help).
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1159, in UnitBase._to(self, other)
1157 return self_decomposed.scale / other_decomposed.scale
-> 1159 raise UnitConversionError(f"'{self!r}' is not a scaled version of '{other!r}'")
UnitConversionError: 'Unit("Mpc2")' is not a scaled version of 'Unit("Mpc")'
During handling of the above exception, another exception occurred:
UnitConversionError Traceback (most recent call last)
Cell In[82], line 18
15 url = 'https://gracedb.ligo.org/api/superevents/S190814bv/files/bayestar.multiorder.fits'
16 skymap = read_sky_map(url, moc=True)
---> 18 result = crossmatch(ex_skymap, gal_coords)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/ligo/skymap/postprocess/crossmatch.py:295, in crossmatch(sky_map, coordinates, contours, areas, modes, cosmology)
293 true_dec = coordinates.dec.rad
294 if np.any(coordinates.distance != 1):
--> 295 true_dist = coordinates.distance.to_value(u.Mpc)
296 else:
297 true_dist = None
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:985, in Quantity.to_value(self, unit, equivalencies)
982 scale = self.unit._to(unit)
983 except Exception:
984 # Short-cut failed; try default (maybe equivalencies help).
--> 985 value = self._to_value(unit, equivalencies)
986 else:
987 value = self.view(np.ndarray)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:891, in Quantity._to_value(self, unit, equivalencies)
888 equivalencies = self._equivalencies
889 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
890 # Standard path, let unit to do work.
--> 891 return self.unit.to(
892 unit, self.view(np.ndarray), equivalencies=equivalencies
893 )
895 else:
896 # The .to() method of a simple unit cannot convert a structured
897 # dtype, so we work around it, by recursing.
898 # TODO: deprecate this?
899 # Convert simple to Structured on initialization?
900 result = np.empty_like(self.view(np.ndarray))
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1195, in UnitBase.to(self, other, value, equivalencies)
1193 return UNITY
1194 else:
-> 1195 return self._get_converter(Unit(other), equivalencies)(value)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1124, in UnitBase._get_converter(self, other, equivalencies)
1121 else:
1122 return lambda v: b(converter(v))
-> 1124 raise exc
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1107, in UnitBase._get_converter(self, other, equivalencies)
1105 # if that doesn't work, maybe we can do it with equivalencies?
1106 try:
-> 1107 return self._apply_equivalencies(
1108 self, other, self._normalize_equivalencies(equivalencies)
1109 )
1110 except UnitsError as exc:
1111 # Last hope: maybe other knows how to do it?
1112 # We assume the equivalencies have the unit itself as first item.
1113 # TODO: maybe better for other to have a `_back_converter` method?
1114 if hasattr(other, "equivalencies"):
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1085, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
1082 unit_str = get_err_str(unit)
1083 other_str = get_err_str(other)
-> 1085 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")
UnitConversionError: 'Mpc2' (area) and 'Mpc' (length) are not convertible
```https://git.ligo.org/lscsoft/ligo.skymap/-/issues/42Using ligo.skymap.moc.rasterize ends up with Aborted (core dumped)2024-01-03T19:15:35ZVASILEIOS SKLIRISUsing ligo.skymap.moc.rasterize ends up with Aborted (core dumped)Hi,
I am a relatively new user of ligo.skymaps and I am trying to take transform a moc skyampap into a nested probability map of nside 64 (order 6). I try to use rectorize function for this in the script below. If I use order more than ...Hi,
I am a relatively new user of ligo.skymaps and I am trying to take transform a moc skyampap into a nested probability map of nside 64 (order 6). I try to use rectorize function for this in the script below. If I use order more than 2 the it crashes with:
munmap_chunk(): invalid pointer
Aborted (core dumped)
```python
import numpy as np
from ligo.skymap.moc import rasterize
from astropy.table import Table
skymap = Table.read('https://gracedb.ligo.org/api/superevents/S231226av/files/Bilby.multiorder.fits')
moc_data = np.array([np.array(skymap['UNIQ']),np.array(skymap['PROBDENSITY'])], dtype=[('UNIQ', int), ('PROBDENSITY', float)])
P = rasterize(moc_data, order=2)
```https://git.ligo.org/lscsoft/ligo.skymap/-/issues/40Calculating marganlized distance moments can give incorrect result2023-09-21T20:37:49ZBrandon PiotrzkowskiCalculating marganlized distance moments can give incorrect resultCalculating the margnaticalized distance moments for MS230627t gives results that are inconsistent with values in the header and appear incorrect. File here:
https://gracedb-playground.ligo.org/api/superevents/MS230627t/files/bayestar.m...Calculating the margnaticalized distance moments for MS230627t gives results that are inconsistent with values in the header and appear incorrect. File here:
https://gracedb-playground.ligo.org/api/superevents/MS230627t/files/bayestar.multiorder.fits,1
```
from ligo.skymap.io import read_sky_map
import ligo.skymap.distance
skymap_multi = read_sky_map('bayestar-MS230627t.multiorder.fits,1', moc=True)
dist_mean, dist_sigma = ligo.skymap.distance.parameters_to_marginal_moments(skymap_multi['PROBDENSITY'], skymap_multi['DISTMU'], skymap_multi['DISTSIGMA'])
print(dist_mean, dist_sigma)
```
```
7862770.23043921 nan
```
These disagrees with the distance moments in the header:
```
print(skymap_multi.meta['distmean'], skymap_multi.meta['diststd'])
```
```
162.5821854394626 43.94262766367867
```Leo P. SingerLeo P. Singerhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/37Refactor Python implementation of BAYESTAR adaptive grid as a reusable function2023-07-06T20:38:38ZLeo P. SingerRefactor Python implementation of BAYESTAR adaptive grid as a reusable functionThere is a Python implementation of the BAYESTAR adaptive grid algorithm in https://git.ligo.org/lscsoft/ligo.skymap/-/blob/v1.0.7/ligo/skymap/kde.py#L325-348.
Refactor it as a standalone function that takes as arguments:
- the initial ...There is a Python implementation of the BAYESTAR adaptive grid algorithm in https://git.ligo.org/lscsoft/ligo.skymap/-/blob/v1.0.7/ligo/skymap/kde.py#L325-348.
Refactor it as a standalone function that takes as arguments:
- the initial nside value
- the number of rounds
- a callable that returns the (possible vector valued) pixel value for a given coordinate
- a callable that returns the sort key for a pixel value
and returns a multi-order map as an Astropy table.
The function might live in https://git.ligo.org/lscsoft/ligo.skymap/-/blob/main/ligo/skymap/moc.py.Brandon PiotrzkowskiBrandon Piotrzkowskihttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/38Vectorize find_ellipse() function over credible levels2023-06-15T13:30:41ZLeo P. SingerVectorize find_ellipse() function over credible levelsTurns out that this should be pretty straightforward. @brandon.piotrzkowski, you'll find this handy. I'll make it happen.Turns out that this should be pretty straightforward. @brandon.piotrzkowski, you'll find this handy. I'll make it happen.Leo P. SingerLeo P. Singerhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/36Crossmatch tests fail with scipy 1.10.0 due to `multivariate_normal_frozen` b...2023-02-22T00:49:09ZJacopo TissinoCrossmatch tests fail with scipy 1.10.0 due to `multivariate_normal_frozen` behavior changeTests in `ligo/skymap/postprocess/tests/test_crossmatch.py` fail when run with the latest scipy 1.10.0
with the error
```
skymap = cartesian_gaussian_to_skymap(
> 6, cartesian_gaussian.mean, cartesian_gaussian.cov)
E ...Tests in `ligo/skymap/postprocess/tests/test_crossmatch.py` fail when run with the latest scipy 1.10.0
with the error
```
skymap = cartesian_gaussian_to_skymap(
> 6, cartesian_gaussian.mean, cartesian_gaussian.cov)
E AttributeError: 'multivariate_normal_frozen' object has no attribute 'cov'
```
This is due to an [undocumented behavior change](https://github.com/scipy/scipy/issues/17896) in `multivariate_normal_frozen` object in scipy (or better put, on the tests here relying on the `rv.cov` attribute being available, which is not specified in the API).
Fixing this on the scipy side seems easy and might even be incorporated in a new patch version
like 1.10.1, but this depends on scipy so it cannot be guaranteed.
Meanwhile, I would simply advise to **not update scipy to 1.10.0** for the moment.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/35Inconsistent tuple lengths in find_ellipse2022-10-07T21:00:34ZBrandon PiotrzkowskiInconsistent tuple lengths in find_ellipseWhen testing in gwcelery, an error was found in `find_ellipse` where the expected values to unpack were different from expected:
![Screen_Shot_2022-10-06_at_10.35.13_AM](/uploads/9916a6cd54cd7e281c8c2ba05e70c84a/Screen_Shot_2022-10-06_at...When testing in gwcelery, an error was found in `find_ellipse` where the expected values to unpack were different from expected:
![Screen_Shot_2022-10-06_at_10.35.13_AM](/uploads/9916a6cd54cd7e281c8c2ba05e70c84a/Screen_Shot_2022-10-06_at_10.35.13_AM.png)
This appears to occur due to the tuples returned being different lengths, with the standard output being six values
(https://git.ligo.org/lscsoft/ligo.skymap/-/blob/ba63039d857d5097417e5d967d7d6db185327d5f/ligo/skymap/postprocess/ellipse.py#L387)
and the error output being five (https://git.ligo.org/lscsoft/ligo.skymap/-/blob/ba63039d857d5097417e5d967d7d6db185327d5f/ligo/skymap/postprocess/ellipse.py#L371)
This could be fixed by changing the latter output to six np.nan values.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/34bayestar-inject fails with AttributeError with Astropy 5.12022-05-27T18:30:25ZDuncan Macleodduncan.macleod@ligo.orgbayestar-inject fails with AttributeError with Astropy 5.1`bayestar-inject` is incompatible with Astropy 5.1, creating an `AttributeError` when creating the command-line parser:
```pytb
$ python -c "from ligo.skymap.tool.bayestar_inject import main; main()" --help
Traceback (most recent call l...`bayestar-inject` is incompatible with Astropy 5.1, creating an `AttributeError` when creating the command-line parser:
```pytb
$ python -c "from ligo.skymap.tool.bayestar_inject import main; main()" --help
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "/home/duncan/opt/mambaforge/envs/py310/lib/python3.10/site-packages/ligo/skymap/tool/bayestar_inject.py", line 272, in main
p = parser()
File "/home/duncan/opt/mambaforge/envs/py310/lib/python3.10/site-packages/ligo/skymap/tool/bayestar_inject.py", line 232, in parser
'--cosmology', choices=cosmology.parameters.available,
File "/home/duncan/opt/mambaforge/envs/py310/lib/python3.10/site-packages/astropy/cosmology/__init__.py", line 38, in __getattr__
raise AttributeError(f"module {__name__!r} has no attribute {name!r}.")
AttributeError: module 'astropy.cosmology' has no attribute 'parameters'.. Did you mean: 'Parameter'?
```https://git.ligo.org/lscsoft/ligo.skymap/-/issues/33bayestar error trapping2022-03-18T16:24:46ZAlexander Pacebayestar error trappingI had this event page brought to my attention this morning: https://gracedb.ligo.org/events/G368545/view/ which is the preferred event for [S200316b](https://gracedb.ligo.org/superevents/S200316bj/view/), a public superevent.
The last ...I had this event page brought to my attention this morning: https://gracedb.ligo.org/events/G368545/view/ which is the preferred event for [S200316b](https://gracedb.ligo.org/superevents/S200316bj/view/), a public superevent.
The last log message (#59) is traceback from an exception in `bayestar_localize_lvalert.py`. I don't necessarily want to delete the log message itself because of the unofficial policy of maintaining data integrity in GraceDB's record, and because it's not on a public-facing page.
Consider improving error trapping in `ligo.skymap` before they make it into the input of `gracedb-client`.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/29Python 3.10 compatibility will (probably) require building against numpy 1.212022-01-18T14:13:47ZDuncan Macleodduncan.macleod@ligo.orgPython 3.10 compatibility will (probably) require building against numpy 1.21I am presuming that the NumPy team will only build Python 3.10 wheels for the current version, so supporting Python 3.10 in _this_ project will likely require updating the build requirement away from numpy 1.19.3.I am presuming that the NumPy team will only build Python 3.10 wheels for the current version, so supporting Python 3.10 in _this_ project will likely require updating the build requirement away from numpy 1.19.3.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/31Failures when running with numpy 1.22.02022-01-10T18:57:36ZDuncan Macleodduncan.macleod@ligo.orgFailures when running with numpy 1.22.0I am seeing a few failures when running against numpy 1.22.0. To demonstrate I took the liberty of running a pipeline on the main branch, see https://git.ligo.org/lscsoft/ligo.skymap/-/pipelines/338728 for examples of errors like this:
...I am seeing a few failures when running against numpy 1.22.0. To demonstrate I took the liberty of running a pipeline on the main branch, see https://git.ligo.org/lscsoft/ligo.skymap/-/pipelines/338728 for examples of errors like this:
```
Warning, treated as error:
/tmp/tmp.QeXs3Yt3M5/docs/tool/ligo_skymap_plot.rst:14:Exception occurred in plotting ligo_skymap_plot-1
from /tmp/tmp.QeXs3Yt3M5/docs/tool/ligo_skymap_plot.rst:
Traceback (most recent call last):
File "/tmp/tmp.QeXs3Yt3M5/.tox/build_docs/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 57, in _wrapfunc
return bound(*args, **kwds)
TypeError: the resolved dtypes are not compatible with add.accumulate. Resolved (dtype('float32'), dtype('float32'), dtype('float32'))
```https://git.ligo.org/lscsoft/ligo.skymap/-/issues/30Changes to python-ligo-lw will break skymap/util/ilwd.py2021-12-07T19:49:30ZIan HarryChanges to python-ligo-lw will break skymap/util/ilwd.pyHi Leo,
Apologies if you are already aware of this. We noticed that forthcoming changes to python-ligo-lw (https://git.ligo.org/kipp.cannon/python-ligo-lw/-/commit/a292464f1c4742548432b2b0ac635dcc86ac82e1) will break the reading code de...Hi Leo,
Apologies if you are already aware of this. We noticed that forthcoming changes to python-ligo-lw (https://git.ligo.org/kipp.cannon/python-ligo-lw/-/commit/a292464f1c4742548432b2b0ac635dcc86ac82e1) will break the reading code defined in skymap/util/ilwd.py as `IDTypes` will be removed. Not sure if it's easily possible to patch this to work with the changes (I ping @tito-canton who has looked at this on our end as well).
Thanks!
Ianhttps://git.ligo.org/lscsoft/ligo.skymap/-/issues/28bayestar-localize-coincs fails if 'file' utility not installed2021-12-01T02:33:21ZDuncan Macleodduncan.macleod@ligo.orgbayestar-localize-coincs fails if 'file' utility not installedThe ligo.skymap test suite fails when calling `bayestar-localize-coincs` if the `file` (`/usr/bin/file` from the `file` package on Debian) utility isn't installed. This isn't a standard package in `ubuntu:focal` containers and similar.
...The ligo.skymap test suite fails when calling `bayestar-localize-coincs` if the `file` (`/usr/bin/file` from the `file` package on Debian) utility isn't installed. This isn't a standard package in `ubuntu:focal` containers and similar.
The utility is required to work out the filetype for some files, is is possible to just read the magic characters from the file headers instead?https://git.ligo.org/lscsoft/ligo.skymap/-/issues/26Error while writing skymaps with version 0.5.22021-04-08T18:29:28ZTito Dal CantonError while writing skymaps with version 0.5.2```
% bayestar-localize-coincs --f-low 20 coinc-1272790100.1640625-hfixfk.xml.gz coinc-1272790100.1640625-hfixfk.xml.gz
2021-04-06 13:13:40,823 INFO Using 8 OpenMP thread(s)
2021-04-06 13:13:40,823 INFO coinc-1272790100.1640625-hfixfk.xm...```
% bayestar-localize-coincs --f-low 20 coinc-1272790100.1640625-hfixfk.xml.gz coinc-1272790100.1640625-hfixfk.xml.gz
2021-04-06 13:13:40,823 INFO Using 8 OpenMP thread(s)
2021-04-06 13:13:40,823 INFO coinc-1272790100.1640625-hfixfk.xml.gz,coinc-1272790100.1640625-hfixfk.xml.gz:reading input files
2021-04-06 13:13:40,946 WARNING Time slide record is missing for 0, guessing that this is zero-lag
2021-04-06 13:13:40,946 INFO 0:computing sky map
2021-04-06 13:13:41,042 WARNING Truncating PSD at 0.95 of maximum frequency to suppress rolloff artifacts. This option may be removed in the future.
2021-04-06 13:13:41,043 WARNING Truncating PSD at 0.95 of maximum frequency to suppress rolloff artifacts. This option may be removed in the future.
2021-04-06 13:13:41,044 WARNING Truncating PSD at 0.95 of maximum frequency to suppress rolloff artifacts. This option may be removed in the future.
2021-04-06 13:13:41,045 WARNING Template is unspecified; using ER10/O2 uberbank criterion
2021-04-06 13:13:41,045 INFO Selected template: SEOBNRv4_ROM
2021-04-06 13:13:41,231 WARNING Assuming PSD is infinite at 0.015625 Hz because PSD is only sampled down to 18 Hz
2021-04-06 13:13:41,231 WARNING Assuming PSD is infinite at 4095.98 Hz because PSD is only sampled up to 972.75 Hz
2021-04-06 13:13:41,237 WARNING Assuming PSD is infinite at 0.015625 Hz because PSD is only sampled down to 18 Hz
2021-04-06 13:13:41,237 WARNING Assuming PSD is infinite at 4095.98 Hz because PSD is only sampled up to 972.75 Hz
2021-04-06 13:13:41,243 WARNING Assuming PSD is infinite at 0.015625 Hz because PSD is only sampled down to 18 Hz
2021-04-06 13:13:41,243 WARNING Assuming PSD is infinite at 4095.98 Hz because PSD is only sampled up to 972.75 Hz
2021-04-06 13:13:55,109 INFO finished computationally-intensive section in real=14.163s, user=100.320s, sys=0.444s
2021-04-06 13:13:55,109 INFO 0:saving sky map
Traceback (most recent call last):
File "/Users/tito/anaconda3/bin/bayestar-localize-coincs", line 8, in <module>
sys.exit(main())
File "/Users/tito/anaconda3/lib/python3.7/site-packages/ligo/skymap/tool/bayestar_localize_coincs.py", line 187, in main
os.path.join(opts.output, filename), sky_map, nest=True)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/ligo/skymap/io/fits.py", line 384, in write_sky_map
hdulist.writeto(filename, overwrite=True)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/utils/decorators.py", line 535, in wrapper
return function(*args, **kwargs)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py", line 919, in writeto
self.verify(option=output_verify)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/verify.py", line 73, in verify
errs = self._verify(opt)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py", line 1245, in _verify
result = hdu._verify(option)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/hdu/table.py", line 535, in _verify
errs = super()._verify(option=option)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py", line 1587, in _verify
errs = super()._verify(option=option)
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py", line 1074, in _verify
errs.append(card._verify(option))
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/card.py", line 1118, in _verify
keyword, valuecomment = self._split()
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/card.py", line 815, in _split
for card in self._itersubcards():
File "/Users/tito/anaconda3/lib/python3.7/site-packages/astropy/io/fits/card.py", line 1171, in _itersubcards
'Long card images must have CONTINUE cards after '
astropy.io.fits.verify.VerifyError: Long card images must have CONTINUE cards after the first card.
```
I am using Python 3.7.7, astropy 4.2.1, ligo-skymap 0.5.2 and numpy 1.20.2. The same command completes successfully with astropy 4.2, ligo.skymap 0.5.0 and numpy 1.19.2.
First reported by Alex Nitz.https://git.ligo.org/lscsoft/ligo.skymap/-/issues/25Incompatibility with healpy in 0.5.02021-02-16T13:38:07ZTito Dal CantonIncompatibility with healpy in 0.5.0Trying to read a FITS skymap produced by ligo.skymap 0.5.0 with healpy leads to the following error.
```
In [1]: import healpy as hp
In [2]:...Trying to read a FITS skymap produced by ligo.skymap 0.5.0 with healpy leads to the following error.
```
In [1]: import healpy as hp
In [2]: hp.__version__
Out[2]: '1.14.0'
In [3]: x = hp.read_map('test.fits.gz')
/Users/tito/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:369: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
"If you are not specifying the input dtype and using the default "
/Users/tito/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 256
warnings.warn("NSIDE = {0:d}".format(nside))
/Users/tito/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = NESTED in fits file
warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-3-81c33a5c1363> in <module>
----> 1 x = hp.read_map('test.fits.gz')
~/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py in read_map(filename, field, dtype, nest, partial, hdu, h, verbose, memmap)
404
405 # partial sky: check OBJECT, then INDXSCHM
--> 406 obj = fits_hdu.header.get("OBJECT", "UNDEF").strip()
407 if obj != "UNDEF":
408 if obj == "PARTIAL":
AttributeError: 'int' object has no attribute 'strip'
```
The issue seems to be that ligo.skymap sets the `OBJECT` header field to an integer, while healpy assumes it is a string-like type.
Would it be a problem to simply convert the integer to a string?https://git.ligo.org/lscsoft/ligo.skymap/-/issues/24Shifted SNR distribution when realizing an event with Gaussian noise multiple...2021-02-07T20:17:20ZGeoffrey MoShifted SNR distribution when realizing an event with Gaussian noise multiple timescc: @sylvia.biscoveanu
When running a matched-filter search on an event multiple times with Gaussian noise, we expect to see each interferometer's SNR distributed as a Gaussian about the optimal SNR with unit variance (see fig. 3 in [1]...cc: @sylvia.biscoveanu
When running a matched-filter search on an event multiple times with Gaussian noise, we expect to see each interferometer's SNR distributed as a Gaussian about the optimal SNR with unit variance (see fig. 3 in [1]).
For a network of N (in this case, four) detectors, the resulting network SNR^2 should be distributed as a noncentral chi-square distribution with N degrees of freedom and the non-centrality parameter equal to the network optimal SNR^2. (see [2])
However, this isn't the distribution that I'm seeing, for 100 matched-filter searches of the same event using `bayestar-realize-coincs`:
![image](/uploads/fbc6c067814f4b5b344a21f3956316c7/image.png)
Here is code to replicate this example. The exact files can also be found on the CIT cluster at `/home/geoffrey.mo/bayestar-test/noise_snr`. My ligo.skymap version is 0.1.16.
```
lalapps_inspinj -o inj.xml --m-distr fixMasses --fixed-mass1 1.3 --fixed-mass2 1.3 --t-distr uniform --time-step 100 --gps-start-time 1000000000 --gps-end-time 1000000099 --d-distr volume --min-distance 150000 --max-distance 160000 --l-distr random --i-distr uniform --f-lower 20 --disable-spin --waveform IMRPhenom_NRTidal
# Make a copy of this, call it inj_many.xml.
cp inj.xml inj_many.xml
# DO THIS MANUALLY #
# Then copy the line in the sim_inspiral inj_many.xml table many times (in this case, 100).
# Make sure to add the delimiter at the end of the line before copying.
bayestar-sample-model-psd -o psd.xml --H1=aLIGOAdVO4T1800545 --L1=aLIGOAdVO4T1800545 --V1=AdVMidLowSensitivityP1200087 --K1=KAGRAEarlySensitivityT1600593
# Get the zero-noise (optimal) SNR:
bayestar-realize-coincs -o coinc.xml inj.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64
# Now localize the many injs:
bayestar-realize-coincs -o coinc_many.xml inj_many.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64 --measurement-error gaussian-noise
# Run this python script to parse the xmls for the SNRs and plot the distribution
python parse_plot.py
```
And this is parse_plot.py:
```
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import xml.etree.ElementTree as ET
# These don't seem to work due to some ligolw incompatabilities with the
# coinc_event_id column validity. I really don't want to debug this right now
# so I'm just going to do some ugly parsing of the xml to get the SNRs.
# coinc_one = EventTable.read('coinc.xml', tablename='coinc_inspiral')
# coinc_many = EventTable.read('coinc_many.xml', tablename='coinc_inspiral')
single_tree = ET.parse('coinc.xml')
single_root = single_tree.getroot()
optimal_snr = float(single_root[1][-1].text.split(',')[-1].split('\n')[0])
many_tree = ET.parse('coinc_many.xml')
many_root = many_tree.getroot()
events = [event.split(',') for event in many_root[1][-1].text.split('\n')[:-1] if
len(event) > 1]
events[-1].append('') # last row needs fixing to look like the other rows
snrs = np.array([float(event[-2]) for event in events])
xarray = np.linspace(5, 15, 200)
plt.figure(figsize=(12,8))
# this is the distribution we expect for the quadrature sum
# of four Gaussians (one for each detector)
plt.hist(np.random.noncentral_chisquare(4, optimal_snr**2, 100000),
bins=200, histtype='step', density=True, label='Noncentral chi squared')
plt.hist(snrs**2, bins=20, histtype='step', density=True, label='Observed');
plt.axvline(optimal_snr**2, ls='--', color='black', label=r'$\rho_\mathrm{opt}$');
plt.legend(loc='upper left');
plt.xlabel('SNR$^2$');
plt.ylabel('Probability density');
plt.savefig('snr2_dist.png', dpi=300)
```
[1] https://arxiv.org/pdf/1809.02293.pdf
[2] Mathematica code:
```
In[1] = TransformedDistribution[
s^2 + t^2 + u^2 + v^2, {s \[Distributed] NormalDistribution[w, 1],
t \[Distributed] NormalDistribution[x, 1],
u \[Distributed] NormalDistribution[y, 1],
v \[Distributed] NormalDistribution[z, 1]}]
Out[1] = NoncentralChiSquareDistribution[4, w^2 + x^2 + y^2 + z^2]
```