Skip to content

Commit

Permalink
Merge eb4a8c7 into 3f10025
Browse files Browse the repository at this point in the history
  • Loading branch information
segasai committed Feb 18, 2024
2 parents 3f10025 + eb4a8c7 commit c743f43
Show file tree
Hide file tree
Showing 11 changed files with 512 additions and 1,226 deletions.
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
83 changes: 39 additions & 44 deletions py/dynesty/bounding.py
Original file line number Diff line number Diff line change
Expand Up @@ -661,43 +661,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 +711,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 +720,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 +745,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 +759,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 +853,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 @@ -927,44 +924,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 +976,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 +1011,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 +1025,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 +1117,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

0 comments on commit c743f43

Please sign in to comment.