Skip to content

Commit

Permalink
Add a Fisher distribution for sky locations (gwastro#4209)
Browse files Browse the repository at this point in the history
* Update transforms.py

Adding new_z_to_euler and rotate_euler from pylal.sphericalutils

* Update sky_location.py

Added new class FisherDist for drawing from sky locations (in ra, dec)

* Update sky_location.py

Included help message for fisherDist

* Update sky_location.py

* Code climate issues

* Code climate issues

* Code-climate changes

* Code climate issues

* Code climate issues

* Update transforms.py

changed lal.PI to numpy.pi in new_z_to_euler

* Update transforms.py

minor changes

* Update sky_location.py

updated docstring, did changes to make the code look more simple

* Update sky_location.py

typos corrected

* Update sky_location.py

code climate changes

* Update sky_location.py

* Update sky_location.py

Included the definitions rvs_polaz and rvs_radec

* Update transforms.py

Included the definitions decra2polaz and polaz2decra

* Update sky_location.py

* Update transforms.py

code climate issues fixed

* Update sky_location.py

code climate issues fixed

* Update sky_location.py

code climate issues

* Update sky_location.py

fixed code climate issues

* Update sky_location.py

fixed code climate issues

* Update sky_location.py

* Update pycbc/distributions/sky_location.py

Co-authored-by: Tito Dal Canton <[email protected]>

* Update sky_location.py

* Update sky_location.py

* Update transforms.py

modified polaz2decra to polaz2radec

* Update sky_location.py

* Update sky_location.py

Forgotten polaz2decra --> polaz2radec

* Update sky_location.py

Code climate fixes

* Update sky_location.py

Updated the code for checking the effectiveness of running from config file

* Create sky_location_old.py

* modified to include from config

* Modified to run from config file

* Added Fisher distribution to the lista of available distributions

* Delete sky_location_old.py

* Update sky_location.py

* Renamed Fisher to FisherSky

* Updated the name to FisherSky

* changed Fisher to FisherSky in __all__

* removed decra2polaz and polaz2radec to conversions

* Included def decra2polaz, polaz2radec in line 859

* Added angle_unit option

* removed decra2polaz and polaz2decra

* made carresction

* made corrections

* Corrected the typing error while removing recra2polaz

* Modified angle_unit and updated init file

* Empty line on 852 retrieved

* corrected

* Empty line on852 retrieved

* Fixed code climate issues

* Sigma conversion taken care of

* Fixed code climate issues

* Update sky_location.py

Codeclimate

* Update test_distributions.py

Adding fisher_sky to the list of distributions that do not undergo the distribution unittests.

* Update sky_location.py

Syntax fixed.

* Update sky_location.py complying to codeclimate

* Refactor using scipy's Rotation

* Improve docstring and make warning stricter

* Update pycbc/distributions/sky_location.py

Co-authored-by: Francesco Pannarale <[email protected]>

* Update pycbc/distributions/sky_location.py

right angle to correct angle

Co-authored-by: Francesco Pannarale <[email protected]>

* Update sky_location.py

logging issue solved

---------

Co-authored-by: Francesco Pannarale <[email protected]>
Co-authored-by: Tito Dal Canton <[email protected]>
Co-authored-by: P Prasia <[email protected]>
Co-authored-by: Tito Dal Canton <[email protected]>
  • Loading branch information
5 people authored and PRAVEEN-mnl committed Jun 19, 2023
1 parent 34d6a21 commit f49620a
Show file tree
Hide file tree
Showing 4 changed files with 147 additions and 6 deletions.
2 changes: 1 addition & 1 deletion pycbc/conversions.py
Original file line number Diff line number Diff line change
Expand Up @@ -862,7 +862,6 @@ def optimal_ra_from_detector(detector_name, tc):
"""For a given detector and GPS time, return the optimal orientation
(directly overhead of the detector) in right ascension.
Parameters
----------
detector_name : string
Expand All @@ -877,6 +876,7 @@ def optimal_ra_from_detector(detector_name, tc):
"""
return optimal_orientation_from_detector(detector_name, tc)[0]


#
# =============================================================================
#
Expand Down
5 changes: 3 additions & 2 deletions pycbc/distributions/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
from pycbc.distributions.arbitrary import Arbitrary, FromFile
from pycbc.distributions.gaussian import Gaussian
from pycbc.distributions.power_law import UniformPowerLaw, UniformRadius
from pycbc.distributions.sky_location import UniformSky
from pycbc.distributions.sky_location import UniformSky, FisherSky
from pycbc.distributions.uniform import Uniform
from pycbc.distributions.uniform_log import UniformLog10
from pycbc.distributions.spins import IndependentChiPChiEff
Expand Down Expand Up @@ -60,7 +60,8 @@
DistributionFunctionFromFile.name: DistributionFunctionFromFile,
FixedSamples.name: FixedSamples,
MchirpfromUniformMass1Mass2.name: MchirpfromUniformMass1Mass2,
QfromUniformMass1Mass2.name: QfromUniformMass1Mass2
QfromUniformMass1Mass2.name: QfromUniformMass1Mass2,
FisherSky.name: FisherSky
}

def read_distributions_from_config(cp, section="prior"):
Expand Down
145 changes: 142 additions & 3 deletions pycbc/distributions/sky_location.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,18 @@
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

"""This modules provides classes for evaluating sky distributions in
right acension and declination.
right ascension and declination.
"""


import logging
import numpy
from scipy.spatial.transform import Rotation
from pycbc.distributions import angular
from pycbc import VARARGS_DELIM
from pycbc.io import FieldArray


class UniformSky(angular.UniformSolidAngle):
"""A distribution that is uniform on the sky. This is the same as
Expand All @@ -32,4 +38,137 @@ class UniformSky(angular.UniformSolidAngle):
_default_azimuthal_angle = 'ra'


__all__ = ['UniformSky']
class FisherSky():
"""A distribution that returns a random angle drawn from an approximate
`Von_Mises-Fisher distribution`_. Assumes that the Fisher concentration
parameter is large, so that we can draw the samples from a simple
rotationally-invariant distribution centered at the North Pole (which
factors as a uniform distribution for the right ascension, and a Rayleigh
distribution for the declination, as described in
`Fabrycky and Winn 2009 ApJ 696 1230`) and then rotate the samples to be
centered around the specified mean position. As in UniformSky, the
declination varies from π/2 to -π/2 and the right ascension varies from
0 to 2π.
.. _Von_Mises-Fisher distribution:
http://en.wikipedia.org/wiki/Von_Mises-Fisher_distribution
.. _Fabrycky and Winn 2009 ApJ 696 1230:
https://doi.org/10.1088/0004-637X/696/2/1230
.. _Briggs et al 1999 ApJS 122 503:
https://doi.org/10.1086/313221
Parameters
----------
mean_ra: float
RA of the center of the distribution.
mean_dec: float
Declination of the center of the distribution.
sigma: float
Spread of the distribution. For the precise interpretation, see Eq 8
of `Briggs et al 1999 ApJS 122 503`_. This should be smaller than
about 20 deg for the approximation to be valid.
angle_unit: str
Unit for the angle parameters: either "deg" or "rad".
"""
name = 'fisher_sky'
_params = ['ra', 'dec']

def __init__(self, **params):
if params['angle_unit'] not in ['deg', 'rad']:
raise ValueError("Only deg or rad is allowed as angle unit")
mean_ra = params['mean_ra']
mean_dec = params['mean_dec']
sigma = params['sigma']
if params['angle_unit'] == 'deg':
mean_ra = numpy.deg2rad(mean_ra)
mean_dec = numpy.deg2rad(mean_dec)
sigma = numpy.deg2rad(sigma)
if mean_ra < 0 or mean_ra > 2 * numpy.pi:
raise ValueError(
f'The mean RA must be between 0 and 2π, {mean_ra} rad given'
)
if mean_dec < -numpy.pi/2 or mean_dec > numpy.pi/2:
raise ValueError(
'The mean declination must be between '
f'-π/2 and π/2, {mean_dec} rad given'
)
if sigma <= 0 or sigma > 2 * numpy.pi:
raise ValueError(
'Sigma must be positive and smaller than 2π '
'(preferably much smaller)'
)
if sigma > 0.35:
logging.warning(
'Warning: sigma = %s rad is probably too large for the '
'Fisher approximation to be valid', sigma
)
self.rayleigh_scale = 0.66 * sigma
# Prepare a rotation that puts the North Pole at the mean position
self.rotation = Rotation.from_euler(
'yz',
[numpy.pi / 2 - mean_dec, mean_ra]
)

@property
def params(self):
return self._params

@classmethod
def from_config(cls, cp, section, variable_args):
tag = variable_args
variable_args = variable_args.split(VARARGS_DELIM)
if set(variable_args) != set(cls._params):
raise ValueError("Not all parameters used by this distribution "
"included in tag portion of section name")
mean_ra = float(cp.get_opt_tag(section, 'mean_ra', tag))
mean_dec = float(cp.get_opt_tag(section, 'mean_dec', tag))
sigma = float(cp.get_opt_tag(section, 'sigma', tag))
angle_unit = cp.get_opt_tag(section, 'angle_unit', tag)
return cls(
mean_ra=mean_ra,
mean_dec=mean_dec,
sigma=sigma,
angle_unit=angle_unit
)

def rvs(self, size):
# Draw samples from a distribution centered on the North pole
np_ra = numpy.random.uniform(
low=0,
high=(2*numpy.pi),
size=size
)
np_dec = numpy.random.rayleigh(
scale=self.rayleigh_scale,
size=size
)

# Convert the samples to intermediate cartesian representation
np_cart = numpy.empty(shape=(size, 3))
np_cart[:, 0] = numpy.cos(np_ra) * numpy.sin(np_dec)
np_cart[:, 1] = numpy.sin(np_ra) * numpy.sin(np_dec)
np_cart[:, 2] = numpy.cos(np_dec)

# Rotate the samples according to our pre-built rotation
rot_cart = self.rotation.apply(np_cart)

# Convert the samples back to spherical coordinates.
# Some unpleasant conditional operations are needed
# to get the correct angle convention.
rot_radec = FieldArray(
size,
dtype=[
('ra', '<f8'),
('dec', '<f8')
]
)
rot_radec['ra'] = numpy.arctan2(rot_cart[:, 1], rot_cart[:, 0])
neg_mask = rot_radec['ra'] < 0
rot_radec['ra'][neg_mask] += 2 * numpy.pi
rot_radec['dec'] = numpy.arcsin(rot_cart[:, 2])
return rot_radec


__all__ = ['UniformSky', 'FisherSky']
1 change: 1 addition & 0 deletions test/test_distributions.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
# some of these distributons have their own specific unit test
EXCLUDE_DIST_NAMES = ["fromfile", "arbitrary",
"external", "external_func_fromfile",
"fisher_sky",
"uniform_solidangle", "uniform_sky",
"independent_chip_chieff",
"uniform_f0_tau", "fixed_samples"]
Expand Down

0 comments on commit f49620a

Please sign in to comment.