import logging
import numpy
from pycbc.tmpltbank.coord_utils import get_cov_params
logger = logging.getLogger('pycbc.tmpltbank.brute_force_methods')
[docs]
def get_physical_covaried_masses(xis, bestMasses, bestXis, req_match,
massRangeParams, metricParams, fUpper,
giveUpThresh = 5000):
"""
This function takes the position of a point in the xi parameter space and
iteratively finds a close point in the physical coordinate space (masses
and spins).
Parameters
-----------
xis : list or array
Desired position of the point in the xi space. If only N values are
provided and the xi space's dimension is larger then it is assumed that
*any* value in the remaining xi coordinates is acceptable.
bestMasses : list
Contains [totalMass, eta, spin1z, spin2z]. Is a physical position
mapped to xi coordinates in bestXis that is close to the desired point.
This is aimed to give the code a starting point.
bestXis : list
Contains the position of bestMasses in the xi coordinate system.
req_match : float
Desired maximum mismatch between xis and the obtained point. If a point
is found with mismatch < req_match immediately stop and return that
point. A point with this mismatch will not always be found.
massRangeParams : massRangeParameters instance
Instance holding all the details of mass ranges and spin ranges.
metricParams : metricParameters instance
Structure holding all the options for construction of the metric
and the eigenvalues, eigenvectors and covariance matrix
needed to manipulate the space.
fUpper : float
The value of fUpper that was used when obtaining the xi_i
coordinates. This lets us know how to rotate potential physical points
into the correct xi_i space. This must be a key in metricParams.evals,
metricParams.evecs and metricParams.evecsCV
(ie. we must know how to do the transformation for
the given value of fUpper)
giveUpThresh : int, optional (default = 5000)
The program will try this many iterations. If no close matching point
has been found after this it will give up.
Returns
--------
mass1 : float
The heavier mass of the obtained point.
mass2 : float
The smaller mass of the obtained point
spin1z : float
The heavier bodies spin of the obtained point.
spin2z : float
The smaller bodies spin of the obtained point.
count : int
How many iterations it took to find the point. For debugging.
mismatch : float
The mismatch between the obtained point and the input xis.
new_xis : list
The position of the point in the xi space
"""
# TUNABLE PARAMETERS GO HERE!
# This states how far apart to scatter test points in the first proposal
origScaleFactor = 1
# Set up
xi_size = len(xis)
scaleFactor = origScaleFactor
bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
count = 0
unFixedCount = 0
currDist = 100000000000000000
while(1):
# If we are a long way away we use larger jumps
if count:
if currDist > 1 and scaleFactor == origScaleFactor:
scaleFactor = origScaleFactor*10
# Get a set of test points with mass -> xi mappings
totmass, eta, spin1z, spin2z, mass1, mass2, new_xis = \
get_mass_distribution([bestChirpmass, bestMasses[1], bestMasses[2],
bestMasses[3]],
scaleFactor, massRangeParams, metricParams,
fUpper)
cDist = (new_xis[0] - xis[0])**2
for j in range(1,xi_size):
cDist += (new_xis[j] - xis[j])**2
if (cDist.min() < req_match):
idx = cDist.argmin()
scaleFactor = origScaleFactor
new_xis_list = [new_xis[ldx][idx] for ldx in range(len(new_xis))]
return mass1[idx], mass2[idx], spin1z[idx], spin2z[idx], count, \
cDist.min(), new_xis_list
if (cDist.min() < currDist):
idx = cDist.argmin()
bestMasses[0] = totmass[idx]
bestMasses[1] = eta[idx]
bestMasses[2] = spin1z[idx]
bestMasses[3] = spin2z[idx]
bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
currDist = cDist.min()
unFixedCount = 0
scaleFactor = origScaleFactor
count += 1
unFixedCount += 1
if unFixedCount > giveUpThresh:
# Stop at this point
diff = (bestMasses[0]*bestMasses[0] * (1-4*bestMasses[1]))**0.5
mass1 = (bestMasses[0] + diff)/2.
mass2 = (bestMasses[0] - diff)/2.
new_xis_list = [new_xis[ldx][0] for ldx in range(len(new_xis))]
return mass1, mass2, bestMasses[2], bestMasses[3], count, \
currDist, new_xis_list
if not unFixedCount % 100:
scaleFactor *= 2
if scaleFactor > 64:
scaleFactor = 1
# Shouldn't be here!
raise RuntimeError
[docs]
def get_mass_distribution(bestMasses, scaleFactor, massRangeParams,
metricParams, fUpper,
numJumpPoints=100, chirpMassJumpFac=0.0001,
etaJumpFac=0.01, spin1zJumpFac=0.01,
spin2zJumpFac=0.01):
"""
Given a set of masses, this function will create a set of points nearby
in the mass space and map these to the xi space.
Parameters
-----------
bestMasses : list
Contains [ChirpMass, eta, spin1z, spin2z]. Points will be placed around
tjos
scaleFactor : float
This parameter describes the radius away from bestMasses that points
will be placed in.
massRangeParams : massRangeParameters instance
Instance holding all the details of mass ranges and spin ranges.
metricParams : metricParameters instance
Structure holding all the options for construction of the metric
and the eigenvalues, eigenvectors and covariance matrix
needed to manipulate the space.
fUpper : float
The value of fUpper that was used when obtaining the xi_i
coordinates. This lets us know how to rotate potential physical points
into the correct xi_i space. This must be a key in metricParams.evals,
metricParams.evecs and metricParams.evecsCV
(ie. we must know how to do the transformation for
the given value of fUpper)
numJumpPoints : int, optional (default = 100)
The number of points that will be generated every iteration
chirpMassJumpFac : float, optional (default=0.0001)
The jump points will be chosen with fractional variation in chirpMass
up to this multiplied by scaleFactor.
etaJumpFac : float, optional (default=0.01)
The jump points will be chosen with fractional variation in eta
up to this multiplied by scaleFactor.
spin1zJumpFac : float, optional (default=0.01)
The jump points will be chosen with absolute variation in spin1z up to
this multiplied by scaleFactor.
spin2zJumpFac : float, optional (default=0.01)
The jump points will be chosen with absolute variation in spin2z up to
this multiplied by scaleFactor.
Returns
--------
Totmass : numpy.array
Total mass of the resulting points
Eta : numpy.array
Symmetric mass ratio of the resulting points
Spin1z : numpy.array
Spin of the heavier body of the resulting points
Spin2z : numpy.array
Spin of the smaller body of the resulting points
Diff : numpy.array
Mass1 - Mass2 of the resulting points
Mass1 : numpy.array
Mass1 (mass of heavier body) of the resulting points
Mass2 : numpy.array
Mass2 (mass of smaller body) of the resulting points
new_xis : list of numpy.array
Position of points in the xi coordinates
"""
# FIXME: It would be better if rejected values could be drawn from the
# full possible mass/spin distribution. However speed in this function is
# a major factor and must be considered.
bestChirpmass = bestMasses[0]
bestEta = bestMasses[1]
bestSpin1z = bestMasses[2]
bestSpin2z = bestMasses[3]
# Firstly choose a set of values for masses and spins
chirpmass = bestChirpmass * (1 - (numpy.random.random(numJumpPoints)-0.5) \
* chirpMassJumpFac * scaleFactor )
etaRange = massRangeParams.maxEta - massRangeParams.minEta
currJumpFac = etaJumpFac * scaleFactor
if currJumpFac > etaRange:
currJumpFac = etaRange
eta = bestEta * ( 1 - (numpy.random.random(numJumpPoints) - 0.5) \
* currJumpFac)
maxSpinMag = max(massRangeParams.maxNSSpinMag, massRangeParams.maxBHSpinMag)
minSpinMag = min(massRangeParams.maxNSSpinMag, massRangeParams.maxBHSpinMag)
# Note that these two are cranged by spinxzFac, *not* spinxzFac/spinxz
currJumpFac = spin1zJumpFac * scaleFactor
if currJumpFac > maxSpinMag:
currJumpFac = maxSpinMag
# Actually set the new spin trial points
if massRangeParams.nsbhFlag or (maxSpinMag == minSpinMag):
curr_spin_1z_jump_fac = currJumpFac
curr_spin_2z_jump_fac = currJumpFac
# Check spins aren't going to be unphysical
if currJumpFac > massRangeParams.maxBHSpinMag:
curr_spin_1z_jump_fac = massRangeParams.maxBHSpinMag
if currJumpFac > massRangeParams.maxNSSpinMag:
curr_spin_2z_jump_fac = massRangeParams.maxNSSpinMag
spin1z = bestSpin1z + ( (numpy.random.random(numJumpPoints) - 0.5) \
* curr_spin_1z_jump_fac)
spin2z = bestSpin2z + ( (numpy.random.random(numJumpPoints) - 0.5) \
* curr_spin_2z_jump_fac)
else:
# If maxNSSpinMag is very low (0) and maxBHSpinMag is high we can
# find it hard to place any points. So mix these when
# masses are swapping between the NS and BH.
curr_spin_bh_jump_fac = currJumpFac
curr_spin_ns_jump_fac = currJumpFac
# Check spins aren't going to be unphysical
if currJumpFac > massRangeParams.maxBHSpinMag:
curr_spin_bh_jump_fac = massRangeParams.maxBHSpinMag
if currJumpFac > massRangeParams.maxNSSpinMag:
curr_spin_ns_jump_fac = massRangeParams.maxNSSpinMag
spin1z = numpy.zeros(numJumpPoints, dtype=float)
spin2z = numpy.zeros(numJumpPoints, dtype=float)
split_point = int(numJumpPoints/2)
# So set the first half to be at least within the BH range and the
# second half to be at least within the NS range
spin1z[:split_point] = bestSpin1z + \
( (numpy.random.random(split_point) - 0.5)\
* curr_spin_bh_jump_fac)
spin1z[split_point:] = bestSpin1z + \
( (numpy.random.random(numJumpPoints-split_point) - 0.5)\
* curr_spin_ns_jump_fac)
spin2z[:split_point] = bestSpin2z + \
( (numpy.random.random(split_point) - 0.5)\
* curr_spin_bh_jump_fac)
spin2z[split_point:] = bestSpin2z + \
( (numpy.random.random(numJumpPoints-split_point) - 0.5)\
* curr_spin_ns_jump_fac)
# Point[0] is always set to the original point
chirpmass[0] = bestChirpmass
eta[0] = bestEta
spin1z[0] = bestSpin1z
spin2z[0] = bestSpin2z
# Remove points where eta becomes unphysical
eta[eta > massRangeParams.maxEta] = massRangeParams.maxEta
if massRangeParams.minEta:
eta[eta < massRangeParams.minEta] = massRangeParams.minEta
else:
eta[eta < 0.0001] = 0.0001
# Total mass, masses and mass diff
totmass = chirpmass / (eta**(3./5.))
diff = (totmass*totmass * (1-4*eta))**0.5
mass1 = (totmass + diff)/2.
mass2 = (totmass - diff)/2.
# Check the validity of the spin values
# Do the first spin
if maxSpinMag == 0:
# Shortcut if non-spinning
pass
elif massRangeParams.nsbhFlag or (maxSpinMag == minSpinMag):
# Simple case where I don't have to worry about correlation with mass
numploga = abs(spin1z) > massRangeParams.maxBHSpinMag
spin1z[numploga] = 0
else:
# Do have to consider masses
boundary_mass = massRangeParams.ns_bh_boundary_mass
numploga1 = numpy.logical_and(mass1 >= boundary_mass,
abs(spin1z) <= massRangeParams.maxBHSpinMag)
numploga2 = numpy.logical_and(mass1 < boundary_mass,
abs(spin1z) <= massRangeParams.maxNSSpinMag)
numploga = numpy.logical_or(numploga1, numploga2)
numploga = numpy.logical_not(numploga)
spin1z[numploga] = 0
# Same for the second spin
if maxSpinMag == 0:
# Shortcut if non-spinning
pass
elif massRangeParams.nsbhFlag or (maxSpinMag == minSpinMag):
numplogb = abs(spin2z) > massRangeParams.maxNSSpinMag
spin2z[numplogb] = 0
else:
# Do have to consider masses
boundary_mass = massRangeParams.ns_bh_boundary_mass
numplogb1 = numpy.logical_and(mass2 >= boundary_mass,
abs(spin2z) <= massRangeParams.maxBHSpinMag)
numplogb2 = numpy.logical_and(mass2 < boundary_mass,
abs(spin2z) <= massRangeParams.maxNSSpinMag)
numplogb = numpy.logical_or(numplogb1, numplogb2)
numplogb = numpy.logical_not(numplogb)
spin2z[numplogb] = 0
if (maxSpinMag) and (numploga[0] or numplogb[0]):
raise ValueError("Cannot remove the guide point!")
# And remove points where the individual masses are outside of the physical
# range. Or the total masses are.
# These "removed" points will have metric distances that will be much, much
# larger than any thresholds used in the functions in brute_force_utils.py
# and will always be rejected. An unphysical value cannot be used as it
# would result in unphysical metric distances and cause failures.
totmass[mass1 < massRangeParams.minMass1*0.9999] = 0.0001
totmass[mass1 > massRangeParams.maxMass1*1.0001] = 0.0001
totmass[mass2 < massRangeParams.minMass2*0.9999] = 0.0001
totmass[mass2 > massRangeParams.maxMass2*1.0001] = 0.0001
# There is some numerical error which can push this a bit higher. We do
# *not* want to reject the initial guide point. This error comes from
# Masses -> totmass, eta -> masses conversion, we will have points pushing
# onto the boudaries of the space.
totmass[totmass > massRangeParams.maxTotMass*1.0001] = 0.0001
totmass[totmass < massRangeParams.minTotMass*0.9999] = 0.0001
if massRangeParams.max_chirp_mass:
totmass[chirpmass > massRangeParams.max_chirp_mass*1.0001] = 0.0001
if massRangeParams.min_chirp_mass:
totmass[chirpmass < massRangeParams.min_chirp_mass*0.9999] = 0.0001
if totmass[0] < 0.00011:
raise ValueError("Cannot remove the guide point!")
mass1[totmass < 0.00011] = 0.0001
mass2[totmass < 0.00011] = 0.0001
# Then map to xis
new_xis = get_cov_params(mass1, mass2, spin1z, spin2z,
metricParams, fUpper)
return totmass, eta, spin1z, spin2z, mass1, mass2, new_xis
[docs]
def stack_xi_direction_brute(xis, bestMasses, bestXis, direction_num,
req_match, massRangeParams, metricParams, fUpper,
scaleFactor=0.8, numIterations=3000):
"""
This function is used to assess the depth of the xi_space in a specified
dimension at a specified point in the higher dimensions. It does this by
iteratively throwing points at the space to find maxima and minima.
Parameters
-----------
xis : list or array
Position in the xi space at which to assess the depth. This can be only
a subset of the higher dimensions than that being sampled.
bestMasses : list
Contains [totalMass, eta, spin1z, spin2z]. Is a physical position
mapped to xi coordinates in bestXis that is close to the xis point.
This is aimed to give the code a starting point.
bestXis : list
Contains the position of bestMasses in the xi coordinate system.
direction_num : int
The dimension that you want to assess the depth of (0 = 1, 1 = 2 ...)
req_match : float
When considering points to assess the depth with, only consider points
with a mismatch that is smaller than this with xis.
massRangeParams : massRangeParameters instance
Instance holding all the details of mass ranges and spin ranges.
metricParams : metricParameters instance
Structure holding all the options for construction of the metric
and the eigenvalues, eigenvectors and covariance matrix
needed to manipulate the space.
fUpper : float
The value of fUpper that was used when obtaining the xi_i
coordinates. This lets us know how to rotate potential physical points
into the correct xi_i space. This must be a key in metricParams.evals,
metricParams.evecs and metricParams.evecsCV
(ie. we must know how to do the transformation for
the given value of fUpper)
scaleFactor : float, optional (default = 0.8)
The value of the scale factor that is used when calling
pycbc.tmpltbank.get_mass_distribution.
numIterations : int, optional (default = 3000)
The number of times to make calls to get_mass_distribution when
assessing the maximum/minimum of this parameter space. Making this
smaller makes the code faster, but at the cost of accuracy.
Returns
--------
xi_min : float
The minimal value of the specified dimension at the specified point in
parameter space.
xi_max : float
The maximal value of the specified dimension at the specified point in
parameter space.
"""
# Find minimum
ximin = find_xi_extrema_brute(xis, bestMasses, bestXis, direction_num, \
req_match, massRangeParams, metricParams, \
fUpper, find_minimum=True, \
scaleFactor=scaleFactor, \
numIterations=numIterations)
# Find maximum
ximax = find_xi_extrema_brute(xis, bestMasses, bestXis, direction_num, \
req_match, massRangeParams, metricParams, \
fUpper, find_minimum=False, \
scaleFactor=scaleFactor, \
numIterations=numIterations)
return ximin, ximax
[docs]
def find_xi_extrema_brute(xis, bestMasses, bestXis, direction_num, req_match, \
massRangeParams, metricParams, fUpper, \
find_minimum=False, scaleFactor=0.8, \
numIterations=3000):
"""
This function is used to find the largest or smallest value of the xi
space in a specified
dimension at a specified point in the higher dimensions. It does this by
iteratively throwing points at the space to find extrema.
Parameters
-----------
xis : list or array
Position in the xi space at which to assess the depth. This can be only
a subset of the higher dimensions than that being sampled.
bestMasses : list
Contains [totalMass, eta, spin1z, spin2z]. Is a physical position
mapped to xi coordinates in bestXis that is close to the xis point.
This is aimed to give the code a starting point.
bestXis : list
Contains the position of bestMasses in the xi coordinate system.
direction_num : int
The dimension that you want to assess the depth of (0 = 1, 1 = 2 ...)
req_match : float
When considering points to assess the depth with, only consider points
with a mismatch that is smaller than this with xis.
massRangeParams : massRangeParameters instance
Instance holding all the details of mass ranges and spin ranges.
metricParams : metricParameters instance
Structure holding all the options for construction of the metric
and the eigenvalues, eigenvectors and covariance matrix
needed to manipulate the space.
fUpper : float
The value of fUpper that was used when obtaining the xi_i
coordinates. This lets us know how to rotate potential physical points
into the correct xi_i space. This must be a key in metricParams.evals,
metricParams.evecs and metricParams.evecsCV
(ie. we must know how to do the transformation for
the given value of fUpper)
find_minimum : boolean, optional (default = False)
If True, find the minimum value of the xi direction. If False find the
maximum value.
scaleFactor : float, optional (default = 0.8)
The value of the scale factor that is used when calling
pycbc.tmpltbank.get_mass_distribution.
numIterations : int, optional (default = 3000)
The number of times to make calls to get_mass_distribution when
assessing the maximum/minimum of this parameter space. Making this
smaller makes the code faster, but at the cost of accuracy.
Returns
--------
xi_extent : float
The extremal value of the specified dimension at the specified point in
parameter space.
"""
# Setup
xi_size = len(xis)
bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
if find_minimum:
xiextrema = 10000000000
else:
xiextrema = -100000000000
for _ in range(numIterations):
# Evaluate extrema of the xi direction specified
totmass, eta, spin1z, spin2z, _, _, new_xis = \
get_mass_distribution([bestChirpmass,bestMasses[1],bestMasses[2],
bestMasses[3]],
scaleFactor, massRangeParams, metricParams,
fUpper)
cDist = (new_xis[0] - xis[0])**2
for j in range(1, xi_size):
cDist += (new_xis[j] - xis[j])**2
redCDist = cDist[cDist < req_match]
if len(redCDist):
if not find_minimum:
new_xis[direction_num][cDist > req_match] = -10000000
currXiExtrema = (new_xis[direction_num]).max()
idx = (new_xis[direction_num]).argmax()
else:
new_xis[direction_num][cDist > req_match] = 10000000
currXiExtrema = (new_xis[direction_num]).min()
idx = (new_xis[direction_num]).argmin()
if ( ((not find_minimum) and (currXiExtrema > xiextrema)) or \
(find_minimum and (currXiExtrema < xiextrema)) ):
xiextrema = currXiExtrema
bestMasses[0] = totmass[idx]
bestMasses[1] = eta[idx]
bestMasses[2] = spin1z[idx]
bestMasses[3] = spin2z[idx]
bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
return xiextrema