# Source code for pycbc.distributions.spins

# Copyright (C) 2017 Collin Capano, Chris Biwer
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
# Public License for more details.
#
# 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 spin distributions of CBCs.
"""
import logging
import numpy

from pycbc import conversions
from pycbc.distributions.uniform import Uniform
from pycbc.distributions.angular import UniformAngle
from pycbc.distributions.power_law import UniformPowerLaw
from pycbc.distributions.arbitrary import Arbitrary
from pycbc.distributions.bounded import get_param_bounds_from_config, \
VARARGS_DELIM, BoundedDist

logger = logging.getLogger('pycbc.distributions.spins')

[docs]class IndependentChiPChiEff(Arbitrary):
r"""A distribution such that :math:\chi_{\mathrm{eff}} and
:math:\chi_p are uniform and independent of each other.

To ensure constraints are applied correctly, this distribution produces all
three components of both spins as well as the component masses.

Parameters
----------
mass1 : BoundedDist, Bounds, or tuple
The distribution or bounds to use for mass1. Must be either a
BoundedDist giving the distribution on mass1, or bounds (as
either a Bounds instance or a tuple) giving the minimum and maximum
values to use for mass1. If the latter, a Uniform distribution will
be used.
mass2 : BoundedDist, Bounds, or tuple
The distribution or bounds to use for mass2. Syntax is the same as
mass1.
chi_eff : BoundedDist, Bounds, or tuple; optional
The distribution or bounds to use for :math:chi_eff. Syntax is the
same as mass1, except that None may also be passed. In that case,
(-1, 1) will be used for the bounds. Default is None.
chi_a : BoundedDist, Bounds, or tuple; optional
The distribution or bounds to use for :math:chi_a. Syntax is the
same as mass1, except that None may also be passed. In that case,
(-1, 1) will be used for the bounds. Default is None.
xi_bounds : Bounds or tuple, optional
The bounds to use for :math:\xi_1 and :math:\xi_2. Must be
:math:\in (0, 1). If None (the default), will be (0, 1).
nsamples : int, optional
The number of samples to use for the internal kde. The larger the
number of samples, the more accurate the pdf will be, but the longer
it will take to evaluate. Default is 10000.
seed : int, optional
Seed value to use for the number generator for the kde. The current
random state of numpy will be saved prior to setting the seed. After
the samples are generated, the state will be set back to what it was.
If None provided, will use 0.
"""
name = "independent_chip_chieff"
_params = ['mass1', 'mass2', 'xi1', 'xi2', 'chi_eff', 'chi_a',
'phi_a', 'phi_s']

def __init__(self, mass1=None, mass2=None, chi_eff=None, chi_a=None,
xi_bounds=None, nsamples=None, seed=None):

if isinstance(mass1, BoundedDist):
self.mass1_distr = mass1
else:
self.mass1_distr = Uniform(mass1=mass1)
if isinstance(mass2, BoundedDist):
self.mass2_distr = mass2
else:
self.mass2_distr = Uniform(mass2=mass2)
# chi eff
if isinstance(chi_eff, BoundedDist):
self.chieff_distr = chi_eff
else:
if chi_eff is None:
chi_eff = (-1., 1.)
self.chieff_distr = Uniform(chi_eff=chi_eff)
if isinstance(chi_a, BoundedDist):
self.chia_distr = chi_a
else:
if chi_a is None:
chi_a = (-1., 1.)
self.chia_distr = Uniform(chi_a=chi_a)
# xis
if xi_bounds is None:
xi_bounds = (0, 1.)
if (xi_bounds[0] > 1. or xi_bounds[0] < 0.) or (
xi_bounds[1] > 1. or xi_bounds[1] < 0.):
raise ValueError("xi bounds must be in [0, 1)")
self.xi1_distr = UniformPowerLaw(dim=0.5, xi1=xi_bounds)
self.xi2_distr = UniformPowerLaw(dim=0.5, xi2=xi_bounds)
# the angles
self.phia_distr = UniformAngle(phi_a=(0,2))
self.phis_distr = UniformAngle(phi_s=(0,2))
self.distributions = {'mass1': self.mass1_distr,
'mass2': self.mass2_distr,
'xi1': self.xi1_distr,
'xi2': self.xi2_distr,
'chi_eff': self.chieff_distr,
'chi_a': self.chia_distr,
'phi_a': self.phia_distr,
'phi_s': self.phis_distr}
# create random variables for the kde
if nsamples is None:
nsamples = 1e4
# save the current random state
rstate = numpy.random.get_state()
# set the seed
if seed is None:
seed = 0
numpy.random.seed(seed)
rvals = self.rvs(size=int(nsamples))
# reset the random state back to what it was
numpy.random.set_state(rstate)
bounds = dict(b for distr in self.distributions.values()
for b in distr.bounds.items())
super(IndependentChiPChiEff, self).__init__(mass1=rvals['mass1'],
mass2=rvals['mass2'], xi1=rvals['xi1'], xi2=rvals['xi2'],
chi_eff=rvals['chi_eff'], chi_a=rvals['chi_a'],
phi_a=rvals['phi_a'], phi_s=rvals['phi_s'],
bounds=bounds)

def _constraints(self, values):
"""Applies physical constraints to the given parameter values.

Parameters
----------
values : {arr or dict}
A dictionary or structured array giving the values.

Returns
-------
bool
Whether or not the values satisfy physical
"""
mass1, mass2, phi_a, phi_s, chi_eff, chi_a, xi1, xi2, _ = \
conversions.ensurearray(values['mass1'], values['mass2'],
values['phi_a'], values['phi_s'],
values['chi_eff'], values['chi_a'],
values['xi1'], values['xi2'])
s1x = conversions.spin1x_from_xi1_phi_a_phi_s(xi1, phi_a, phi_s)
s2x = conversions.spin2x_from_mass1_mass2_xi2_phi_a_phi_s(mass1, mass2,
xi2, phi_a, phi_s)
s1y = conversions.spin1y_from_xi1_phi_a_phi_s(xi1, phi_a, phi_s)
s2y = conversions.spin2y_from_mass1_mass2_xi2_phi_a_phi_s(mass1, mass2,
xi2, phi_a, phi_s)
s1z = conversions.spin1z_from_mass1_mass2_chi_eff_chi_a(mass1, mass2,
chi_eff, chi_a)
s2z = conversions.spin2z_from_mass1_mass2_chi_eff_chi_a(mass1, mass2,
chi_eff, chi_a)
test = ((s1x**2. + s1y**2. + s1z**2.) < 1.) & \
((s2x**2. + s2y**2. + s2z**2.) < 1.)
return test

def __contains__(self, params):
"""Determines whether the given values are in each parameter's bounds
and satisfy the constraints.
"""
isin = all([params in dist for dist in self.distributions.values()])
if not isin:
return False
# in the individual distributions, apply constrains
return self._constraints(params)

def _draw(self, size=1, **kwargs):
"""Draws random samples without applying physical constrains.
"""
# draw masses
try:
mass1 = kwargs['mass1']
except KeyError:
mass1 = self.mass1_distr.rvs(size=size)['mass1']
try:
mass2 = kwargs['mass2']
except KeyError:
mass2 = self.mass2_distr.rvs(size=size)['mass2']
# draw angles
try:
phi_a = kwargs['phi_a']
except KeyError:
phi_a = self.phia_distr.rvs(size=size)['phi_a']
try:
phi_s = kwargs['phi_s']
except KeyError:
phi_s = self.phis_distr.rvs(size=size)['phi_s']
# draw chi_eff, chi_a
try:
chi_eff = kwargs['chi_eff']
except KeyError:
chi_eff = self.chieff_distr.rvs(size=size)['chi_eff']
try:
chi_a = kwargs['chi_a']
except KeyError:
chi_a = self.chia_distr.rvs(size=size)['chi_a']
# draw xis
try:
xi1 = kwargs['xi1']
except KeyError:
xi1 = self.xi1_distr.rvs(size=size)['xi1']
try:
xi2 = kwargs['xi2']
except KeyError:
xi2 = self.xi2_distr.rvs(size=size)['xi2']
dtype = [(p, float) for p in self.params]
arr = numpy.zeros(size, dtype=dtype)
arr['mass1'] = mass1
arr['mass2'] = mass2
arr['phi_a'] = phi_a
arr['phi_s'] = phi_s
arr['chi_eff'] = chi_eff
arr['chi_a'] = chi_a
arr['xi1'] = xi1
arr['xi2'] = xi2
return arr

[docs]    def apply_boundary_conditions(self, **kwargs):
return kwargs

[docs]    def rvs(self, size=1, **kwargs):
"""Returns random values for all of the parameters.
"""
size = int(size)
dtype = [(p, float) for p in self.params]
arr = numpy.zeros(size, dtype=dtype)
remaining = size
keepidx = 0
while remaining:
draws = self._draw(size=remaining, **kwargs)
remaining = size - keepidx
return arr

[docs]    @classmethod
def from_config(cls, cp, section, variable_args):
"""Returns a distribution based on a configuration file. The parameters
for the distribution are retrieved from the section titled
"[section-variable_args]" in the config file.

Parameters
----------
cp : pycbc.workflow.WorkflowConfigParser
A parsed configuration file that contains the distribution
options.
section : str
Name of the section in the configuration file.
variable_args : str
The names of the parameters for this distribution, separated by
prior.VARARGS_DELIM. These must appear in the "tag" part
of the section header.

Returns
-------
IndependentChiPChiEff
A distribution instance.
"""
tag = variable_args
variable_args = variable_args.split(VARARGS_DELIM)
if not set(variable_args) == set(cls._params):
raise ValueError("Not all parameters used by this distribution "
"included in tag portion of section name")
# get the bounds for the setable parameters
mass1 = get_param_bounds_from_config(cp, section, tag, 'mass1')
mass2 = get_param_bounds_from_config(cp, section, tag, 'mass2')
chi_eff = get_param_bounds_from_config(cp, section, tag, 'chi_eff')
chi_a = get_param_bounds_from_config(cp, section, tag, 'chi_a')
xi_bounds = get_param_bounds_from_config(cp, section, tag, 'xi_bounds')
if cp.has_option('-'.join([section, tag]), 'nsamples'):
nsamples = int(cp.get('-'.join([section, tag]), 'nsamples'))
else:
nsamples = None
return cls(mass1=mass1, mass2=mass2, chi_eff=chi_eff, chi_a=chi_a,
xi_bounds=xi_bounds, nsamples=nsamples)