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

Parallel uniform sampling #471

Open
wants to merge 24 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
9 changes: 0 additions & 9 deletions docs/source/api.rst
Original file line number Diff line number Diff line change
Expand Up @@ -40,15 +40,6 @@ Baseline Sampler
:private-members:
:show-inheritance:

Static Nested Samplers
======================

.. automodule:: dynesty.nestedsamplers
:members:
:undoc-members:
:private-members:
:show-inheritance:

Dynamic Nested Sampler
======================

Expand Down
4 changes: 1 addition & 3 deletions docs/source/references.rst
Original file line number Diff line number Diff line change
Expand Up @@ -116,9 +116,7 @@ Much of the nested sampling error analysis is based on:
*Properties of Nested Sampling.*
Biometrika, 97, 741.

The nested sampling algorithms in
:class:`~dynesty.nestedsamplers.RadFriendsSampler` and
:class:`~dynesty.nestedsamplers.SupFriendsSampler`
The nested sampling algorithms with cubes, balls bounds
are based on:

`Buchner 2016 <http://adsabs.harvard.edu/abs/2014arXiv1407.5459B>`_.
Expand Down
87 changes: 43 additions & 44 deletions py/dynesty/bounding.py
Original file line number Diff line number Diff line change
Expand Up @@ -638,6 +638,8 @@ def update(self,
class RadFriends:
"""
A collection of N-balls of identical size centered on each live point.
Importantly this class requires that the centers of the balls are set
in the .ctrs attribute.

Parameters
----------
Expand All @@ -661,43 +663,43 @@ def __init__(self, ndim, cov=None):

detsign, detln = linalg.slogdet(self.am)
assert detsign > 0
self.logvol_ball = logvol_prefactor(self.n) - 0.5 * detln
self.logvol = logvol_prefactor(self.n) - 0.5 * detln
self.expand = 1.
self.funit = 1

def scale_to_logvol(self, logvol):
"""Scale ball to encompass a target volume."""

f = np.exp((logvol - self.logvol_ball) * (1.0 / self.n))
f = np.exp((logvol - self.logvol) * (1.0 / self.n))
# linear factor
self.cov *= f**2
self.am /= f**2
self.axes *= f
self.axes_inv /= f
self.logvol_ball = logvol
self.logvol = logvol

def within(self, x, ctrs):
def within(self, x):
"""Check which balls `x` falls within."""

# Execute a brute-force search over all balls.
idxs = np.where(
lalg.norm(np.dot(ctrs - x, self.axes_inv), axis=1) <= 1.)[0]
lalg.norm(np.dot(self.ctrs - x, self.axes_inv), axis=1) <= 1.)[0]

return idxs

def overlap(self, x, ctrs):
def overlap(self, x):
"""Check how many balls `x` falls within."""

q = len(self.within(x, ctrs))
q = len(self.within(x))

return q

def contains(self, x, ctrs):
def contains(self, x):
"""Check if the set of balls contains `x`."""

return self.overlap(x, ctrs) > 0
return self.overlap(x) > 0

def sample(self, ctrs, rstate=None, return_q=False):
def sample(self, rstate=None, return_q=False):
"""
Sample a point uniformly distributed within the *union* of balls.

Expand All @@ -711,7 +713,7 @@ def sample(self, ctrs, rstate=None, return_q=False):

"""

nctrs = len(ctrs) # number of balls
nctrs = len(self.ctrs) # number of balls

while True:
ds = randsphere(self.n, rstate=rstate)
Expand All @@ -720,20 +722,20 @@ def sample(self, ctrs, rstate=None, return_q=False):
# If there is only one ball, sample from it.
if nctrs == 1:
q = 1
x = ctrs[0] + dx
x = self.ctrs[0] + dx
else:
# Select a ball at random.
idx = rstate.integers(nctrs)
x = ctrs[idx] + dx
q = self.overlap(x, ctrs)
x = self.ctrs[idx] + dx
q = self.overlap(x)
# random is faster than uniform
if q == 1 or return_q or rstate.random() < (1. / q):
if return_q:
return x, q
else:
return x

def samples(self, nsamples, ctrs, rstate=None):
def samples(self, nsamples, rstate=None):
"""
Draw `nsamples` samples uniformly distributed within the *union* of
balls.
Expand All @@ -745,13 +747,11 @@ def samples(self, nsamples, ctrs, rstate=None):

"""

xs = np.array(
[self.sample(ctrs, rstate=rstate) for i in range(nsamples)])
xs = np.array([self.sample(rstate=rstate) for i in range(nsamples)])

return xs

def monte_carlo_logvol(self,
ctrs,
ndraws=10000,
rstate=None,
return_overlap=True):
Expand All @@ -761,12 +761,11 @@ def monte_carlo_logvol(self,

# Estimate volume using Monte Carlo integration.
samples = ([
self.sample(ctrs, rstate=rstate, return_q=True)
for i in range(ndraws)
self.sample(rstate=rstate, return_q=True) for i in range(ndraws)
])
qs = np.array([_[1] for _ in samples])
qsum = np.sum(1. / qs)
logvol = np.log(1. / ndraws * qsum * len(ctrs)) + self.logvol_ball
logvol = np.log(1. / ndraws * qsum * len(self.ctrs)) + self.logvol

if return_overlap:
# Estimate the fractional amount of overlap with the
Expand Down Expand Up @@ -856,7 +855,7 @@ def update(self,
# Compute volume.
detsign, detln = linalg.slogdet(self.am)
assert detsign > 0
self.logvol_ball = (logvol_prefactor(self.n) - 0.5 * detln)
self.logvol = (logvol_prefactor(self.n) - 0.5 * detln)
self.expand = 1.

# Estimate the volume and fractional overlap with the unit cube
Expand Down Expand Up @@ -904,6 +903,8 @@ def _get_covariance_from_clusters(self, points):
class SupFriends:
"""
A collection of N-cubes of identical size centered on each live point.
Importantly this class requires that the centers of the cubes are set
in the .ctrs attribute.

Parameters
----------
Expand All @@ -927,44 +928,45 @@ def __init__(self, ndim, cov=None):

detsign, detln = linalg.slogdet(self.am)
assert detsign > 0
self.logvol_cube = self.n * np.log(2.) - 0.5 * detln
self.logvol = self.n * np.log(2.) - 0.5 * detln
self.expand = 1.
self.funit = 1

def scale_to_logvol(self, logvol):
"""Scale cube to encompass a target volume."""

f = np.exp((logvol - self.logvol_cube) * (1.0 / self.n))
f = np.exp((logvol - self.logvol) * (1.0 / self.n))
# linear factor
self.cov *= f**2
self.am /= f**2
self.axes *= f
self.axes_inv /= f
self.logvol_cube = logvol
self.logvol = logvol

def within(self, x, ctrs):
def within(self, x):
"""Checks which cubes `x` falls within."""

# Execute a brute-force search over all cubes.
idxs = np.where(
np.max(np.abs(np.dot(ctrs - x, self.axes_inv)), axis=1) <= 1.)[0]
np.max(np.abs(np.dot(self.ctrs -
x, self.axes_inv)), axis=1) <= 1.)[0]

return idxs

def overlap(self, x, ctrs):
def overlap(self, x):
"""Checks how many cubes `x` falls within, skipping the `j`-th
cube."""

q = len(self.within(x, ctrs))
q = len(self.within(x))

return q

def contains(self, x, ctrs):
def contains(self, x):
"""Checks if the set of cubes contains `x`."""

return self.overlap(x, ctrs) > 0
return self.overlap(x) > 0

def sample(self, ctrs, rstate=None, return_q=False):
def sample(self, rstate=None, return_q=False):
"""
Sample a point uniformly distributed within the *union* of cubes.

Expand All @@ -978,30 +980,30 @@ def sample(self, ctrs, rstate=None, return_q=False):

"""

nctrs = len(ctrs) # number of cubes
nctrs = len(self.ctrs) # number of cubes

while True:
ds = rstate.uniform(-1, 1, size=self.n)
dx = np.dot(ds, self.axes)
# If there is only one cube, sample from it.
if nctrs == 1:
x = ctrs[0] + dx
x = self.ctrs[0] + dx
q = 1
else:
# Select a cube at random.
idx = rstate.integers(nctrs)
x = ctrs[idx] + dx
x = self.ctrs[idx] + dx
# Check how many cubes the point lies within, passing over
# the `idx`-th cube `x` was sampled from.
q = self.overlap(x, ctrs)
q = self.overlap(x)
# random() is faster than uniform()
if q == 1 or return_q or rstate.random() < (1. / q):
if return_q:
return x, q
else:
return x

def samples(self, nsamples, ctrs, rstate=None):
def samples(self, nsamples, rstate=None):
"""
Draw `nsamples` samples uniformly distributed within the *union* of
cubes.
Expand All @@ -1013,13 +1015,11 @@ def samples(self, nsamples, ctrs, rstate=None):

"""

xs = np.array(
[self.sample(ctrs, rstate=rstate) for i in range(nsamples)])
xs = np.array([self.sample(rstate=rstate) for i in range(nsamples)])

return xs

def monte_carlo_logvol(self,
ctrs,
ndraws=10000,
rstate=None,
return_overlap=True):
Expand All @@ -1029,11 +1029,10 @@ def monte_carlo_logvol(self,

# Estimate the volume using Monte Carlo integration.
samples = [
self.sample(ctrs, rstate=rstate, return_q=True)
for i in range(ndraws)
self.sample(rstate=rstate, return_q=True) for i in range(ndraws)
]
qsum = sum((1. / q for (x, q) in samples))
logvol = np.log(1. * qsum / ndraws * len(ctrs)) + self.logvol_cube
logvol = np.log(1. * qsum / ndraws * len(self.ctrs)) + self.logvol

if return_overlap:
# Estimate the fractional overlap with the unit cube using
Expand Down Expand Up @@ -1122,7 +1121,7 @@ def update(self,

detsign, detln = linalg.slogdet(self.am)
assert detsign > 0
self.logvol_cube = (self.n * np.log(2.) - 0.5 * detln)
self.logvol = (self.n * np.log(2.) - 0.5 * detln)
self.expand = 1.

# Estimate the volume and fractional overlap with the unit cube
Expand Down
Loading
Loading